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Preface 


This  research  project  is  a  part  of  the  ongoing  effort 
here  at  the  Air  Force  Institute  of  Technology  to  develop 
alternate  methods  for  solving  the  two-dimensional  steady 
state  anisotropic  neutral  particle  transport  equation.  The 
thrust  of  this  research  effort  is  towards  an  accurate  and 
cost-effective  solution  of  the  steady  state  transport  of 
neutrons,  gamma  rays  and  high  energy  x-rays  from  a  low 
altitude  nuclear  burst.  This  problem  which  is  modeled  as 
a  point  source  in  a  two-dimensional  cylindrical  (r,z) 
geometry  with  the  air  ground  interface  included,  is  of 
particular  interest  in  the  areas  of  nuclear  weapons  effects 
and  radiation  physics. 

Presently  the  most  widely  used  computational  methods 
for  solving  the  (air-over-ground)  problem  are  Monte  Carlo 
and  discrete  ordinates.  However,  these  methods  have 
severe  limitations  and  computational  problems.  My  research 
plan  was  to  formulate  and  evaluate  a  solution  technique 
which  did  not  have  these  disadvantages.  A  finite  element 
solution  method  which  is  based  on  a  space-angle  synthesis 
flux  expansion  of  bicubic  splines  and  spherical  harmonics 
was  chosen.  The  merits  of  this  solution  technique  were 
examined  anu  a  comouter  algorithm  for  the  numerical  solu¬ 
tion  of  this  problem  was  developed. 

T  wish  to  a c Know  ledge  and  express  my  appreciation 
for  the  assistance  and  encouragement  which  I  have  received 
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from  the  staff  and  students  of  the  Air  Force  Institute  of 
Technology.  Special  thanks  are  due  to  my  advisor,  Captain 
David  D.  Hardin,  without  whose  direction,  encouragement  and 
many  hours  of  discussion  and  counselling  this  thesis  would 
not  have  been  possible.  I  am  also  grateful  for  the  support, 
advice  and  encouragement  that  was  provided  by  Dr.  J.  Jones 
of  the  Air  Force  Institute  of  Technology  Mathematics 
Department . 

Finally,  I  wish  to  express  my  appreciation  to  my  wife 
and  daughter  for  their  understanding,  patience  and  constant 
support  throughout  this  project.  To  my  wife,  Cynthia,  I 
must  also  express  a  special  thanks  for  her  effort  in  typing 
this  thesis. 
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Abstract 


A  finite  element  space-angle  synthesis 'Solution  of  the 
steady  state  anisotropic  Boltzmann  (transport)  equation  in  a 
two-dimensiona 1  cylindrical  geometry  has  been  developed. 
Starting  from  a  variational  principle  the  Bubnov-Ga 1 erkin 
solution  method  was  applied  to  the  second  order  even  parity 
form  of  the  Boltzmann  equation.  A  trial  function  flux  ex¬ 
pansion  in  bicubic  splines  and  spherical  (surface)  harmonics 
was  used.  A  first  scatter  (collision)  source  and  an  exponen¬ 
tially  varying  atmosphere  was  also  incorporated  into  this 
development . 

Finite  element  space-angle  synthesis  (FESAS)  was 
developed  as  an  alternate  solution  approach  and  an  im¬ 
provement  in  comparison  to  the  methods  of  Monte  Carlo  and 
discrete  ordinates.  FESAS  does  not  have  the  inherent 
characteristics  which  have  produced  the  ray  effect  problem 
in  discrete  ordinates.  Also,  FESAS  may  result  in  lower 
compu ta t iona 1  costs-  than  those  of  Monte  Carlo  and  discrete 
ordinates. 

The  second  order  even  parity  form  of  the  Boltzmann 
equation  was  derived  and  shown  to  be  symmetric,  positive 
definite  and  self-adjoint.  The  equivalence  of  a  varia¬ 
tional  minimization  principle  and  the  Bubnov-Ga  lerkin 
method  of  weighted  residuals  was  established.  The  finite 
•*  1 orient  space  angle  synthesis  a vs  tom  of  equations  was 
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expanded  and  a  numerical  computer  solution  approach  was  im¬ 
plemented.  A  computer  program  was  written  to  solve  for  the 
trial  function  expansion  (mixing)  coefficients,  and  also  to 
compute  the  particle  flux. 


THE  APPLICATION  OF  FINITE  ELEMENTS  AND 
SPACE-ANCLE  SYNTHESIS  TO  THE  ANISOTROPIC 
STEADY  STATE  BOLTZMANN  (TRANSPORT)  EQUATION 


I  Introduction 


Background 

Tho  Air-Over-Ground  Problem.  The  transport  of  neutrons, 
gamma  rays  and  high-energv  x-rays,  away  from  a  low  altitude 
nuclear  explosion  (air-burst),  is  of  special  interest  in 
assessing  the  vulnerability  and  survivability  of  military 
weapon  systems  and  in  making  radiation  exposure  and  dose 
predictions.  This  neutral  particle  transport  problem  in¬ 
creases  in  complexity  because  of  the  exponentially  varying 
air  density  and  the  air-ground  interface.  A  description  of 
neutral  particle  transport  and,  therefore,  the  air-over- 
grcund  problem  is  given  by  the  Boltzmann  transport  equation. 

Numerical  solutions  to  this  problem  already  exist. 

The  main  solution  techniques  are  Monte  Carlo  and  discrete 
ordinates.  However,  discrete  ordinates  and  Monte  Carlo  have 
severe  difficulties  and  disadvantages.  To  perform  an  accurat 
■'  :*  •_  ■  r ;  e  Ui.it!  rec.u :  res  utio  us-’  e:  hours  of  cost  lv 

■.!  '  1  • :  ’  L  >.’*  O  'C '  *  i  tla  L  < .)  S  U  S  l'  j  LOSS  C  v.  QU  ;  HI  ll  O  T  OX(^" 

cui  i  oil  c  i.:no  ti'.an  Mo:;  Co  Carlo,  however ,  it  is  subject:  to  a 
oempuial i onal  difficulty  called  ray  effects  (Ref  1:357). 


Ray  Effects  and  Discrete  Ordinates.  Ray  effects  are  a 
result  of  the  angular  discretization  of  the  particle  flux  in 
the  discrete  ordinate  method.  It  is  not  a  numerical  problem, 
but  originates  in  the  derivation  of  the  discrete  ordinate  Sn 
equations.  In  a  physical  sense  these  equations  only  allow 
source  particles  to  travel  in  specific  directions.  However, 
in  most  practical  problems  these  particles  move  in  all 
directions.  An  in-depth  analysis  of  the  Sn  equations  and  the 
nature  and  reasons  for  ray  effects  can  be  found  elsewhere 
(Ref  1:357). 

Ray  effects  produce  non-physical  distortions  of  the 
angular  flux  in  regions  where  there  are  strong  absorbers, 
localized  sources,  or  high  energy  streaming  particles  (Kef 
2:255-268).  These  distortions  in  the  numerical  formulation 
of  the  discrete  ordinate  method  produce  solutions  which  are 
inaccurate.  The  degree  of  these  inaccuracies  is  dependent 
upon  the  specific  problem  and  the  nature  of  the  absorbing 
media  and  sources.  In  the  air-over-ground  problem  this 
effect  will  bo  significant  because  it  is  essential lv  a 
streaming  particle  problem  with  localized  first  scatter 
sources  . 

A  considerable  amount  of  work  has  already  been  done 

in  an  attemnt  to  eliminate  rav  effects  from  the  S  equations 

n 

and  the  discrete  ordinate  method.  Ray  effects  can  be  initi¬ 
ated  bv  the  use  of  a  fine  angular  mesh  in  the  finite 

differencing  scheme  of  the  S  equations.  However,  this 

n 

approach  increases  the  computat  ional  time  and  the  already 


high  computer  costs.  Other  approaches  involve  a  spherical 
harmonic-like  formulation  (Ref  2:255-268),  piecewise  bilinear 
finite  element  approximations  (Ref  3:205-217),  and  space- 
angle  synthesis  with  specially  tailored  trial  functions  (Ref 
4:322-343)  . 


The  purpose  of  this  research  project  was  to  d  velop  a 
finite  element  solution  to  the  air-over-ground  problem  by 
using  a  space-angle  synthesis  of  bicubic  splines  in  space 
and  spherical  (surface)  harmonics  in  angle.  Specifically, 
a  solution  of  the  monoenergetic  Boltzmann  equation  in  the 
context  of  the  air-over-ground  problem  was  sought.  Working 
from  a  variational  principle  and  using  a  judicious  choice  of 
trial  functions  the  problem  of  ray  effects  may  be  eliminated 
(Ref  3:214).  Also,  this  judicious  trial  function  choice, 
and  a  Bubnov-Galerkin  solution  method  may  be  more  efficient 
and  less  costly  than  Monte  Carlo  or  discrete  ordinates. 

The  steady  state  solution  of  the  Boltzmann  equation 
with  first  scatter  sources,  anisotropic  scattering  and  an 
exponentially  changing  atmosphere  is  desired.  The  problem 
is  formulated  from  a  variational  principle  and  in  a  two- 
dimensional  cylindrical  (r,z)  spatial  geometry  with  the  air- 
over-ground  interface  included.  Fluence  values  as  a 
function  of  two  spatial  (r,z)  and  two  angular  (g,v)  variables 
are  sought.  This  is  a  four  dimensional  problem.  Finally, 
a  numerical  solution  algorithm  and  the  computer  implemen¬ 
tation  of  this  problem  is  required. 
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Assumptions 

There  are  two  basic  assumptions  which  are  made  in  the 
formula; Lon  of  this  problem.  A  time-independent  (steady 
state)  solution  and  axial  symmetry  is  assumed.  Because  of 
the  exponentially  changing  air  density  the  air-over-ground 
problem  is  four  dimensional  with  a  spatially  dependent 
(r,z)  solution.  An  assumption  of  axial  symmetry  is  made 
possible  by  ignoring  the  curvature  of  the  earth.  Within 
the  problem  domain  of  most  practical  problems  the  curvature 
of  the  earth  is  small  and  can  therefore  be  ignored.  Figure 
1  shows  the  spatial  cylindrical  problem  geometry. 

The  flux  from  an  air  burst  is  non-zero  for  a  fraction 
of  a  second  (microseconds).  Therefore,  particle  fluence 
( number /area )  and  not  flux  is  the  more  useful  quantity. 

A  steady  state  formulation  ol  the  air-over-ground  problem 
is  obtained  by  integrating  the  time  dependence  out  of  the 
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Boltzmann  equation.  This  integration  which  is  carried  out 
over  time  limits  when  the  flux  is  zero  produces  a  time 
integrated  or  fluence  equation. 

Development 

In  the  next  chapter  of  this  report  the  problem  equa¬ 
tions  are  presented.  A  finite  element  formulation  of  the 
air-over-ground  problem  is  developed  in  Chapter  III.  The 
Bubnov-Galerkin  method  of  weighted  residuals  is  incorpora¬ 
ted  into  this  development.  In  Chapter  IV  a  space-angle 
synthesis  of  bicubic  splines  and  spherical  harmonics  is 
performed.  An  interpolation  of  the  source  terms  is  also 
outlined  in  Chapter  IV.  A  computer  implementation  of  the 
problem  solution  is  examined  in  Chapter  V.  Finally,  con¬ 
clusions  and  recommendations  are  presented  in  Chapter  VI. 


II  The  Problem  Equation 


The  application  of  finite  elements  and  a  variational 
principle  to  the  air-over-ground  problem  and  the  mono- 
energetic  steady  state  Boltzmann  equation  is  not  a  new 
concept  (Ref  5).  As  in  the  work  by  Wheaton  (Ref  5)  only 
the  monoenergetic  problem  will  be  considered.  It  is 
assumed  that  energy  dependence  can  be  easily  incorporated 
into  this  treatment  by  the  use  of  standard  multigroup 
methods.  The  air-over-ground  problem  which  is  in  effect 
the  steady  state  transport  of  neutral  particles  can  be 
described  by  the  Boltzmann  (transport)  equation  and 
appropriate  boundary  conditions  as  follows: 


a • V  d ( r , fl )  +  ot(r)$(r,£) 


os(r,r.«ir)i(r,JOdfl  ' 

+  S(r,P.) 


(1) 


This  is  the  one  speed  monoenergetic  Boltzmann  equation  in 
general  geometry  whore 


r  =  the  spatial  position  vector, 

Q  =  a  unit  direction  or  velocity  vector, 

V  =  gradient  oDerator, 

4>  =  angular  particle  fluence  in  particles/ 
unit  area /s t erad ian , 

o1-(r)  =  total  macroscopic  interaction  cross 
section  at  spatial  position  r, 


■j1’  (  r  ,  j  •  /)  =  macroscopic  scattering  cross  section. 

1  he  prooability  of  a  particle  at  position 
r  and  di  rec  t  ion  ..  '  scatterin'  into 
direction  a.  It  is  a  tunc i ton  of  the 
scattering  angle  and  not  a  function 

of  the  individual  directions  (isotropic 
mod  in), 
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O'* 


(Ref  6:57).  u  is  the  cosine  of  the  angle  formed  by  the 
z-axis  and  the  particle  velocity  vector  ft .  y  is  the  angle 
between  the  planes  formed  by  the  ?  vector  and  the  z-axis 
and  that  of  the  ft.  vector  and  z-axis. 

The  scattering  properties  of  air  show  a  directional 
dependence  which  is  highly  peaked  in  the  forward  direction 
especially  at  the  high  particle  energies  that  exist  in  the 
air-over-ground  problem.  Because  of  this  the  exterior  boundary 
condition  for  this  problem  will  be  approximated  by  a  vacuum 
boundary  condition: 

•  (?»!)  =  0  for  rs  on  the  boundary  of  the  problem 

domain  and  ft*n  <  0  ( 2 ) 

where  h  is  the  outward  unit  normal  to  the  boundary  surface. 

In  physical  terms  this  is  a  non-reentrant  boundary  condition. 

No  particles  are  allowed  to  reenter  the  region  once  they 
leave  it. 

In  two-d imens iona 1  cylindrical  (r,z)  geometry  there 
is  svnnetrv  in  the  ancle  this  symmetry  can  be  written 

as 

=  r(r,ft2)  r'or  ^p2=^pl"2(ft*api)ft  '•3a’) 

o  r  a  •> 

:(r. C  3  b  > 


3 


4 


This  symmetry  in  >;  is  shown  in  Figure  3,  where  the  vectors 


h  and  6  are  perpendicular. 

Note  that  because  of  the  exponentially  varying  air 
density  (in  the  z  direction),  only  azimuthal  symmetry  in  the 


angle  y  is  assured.  There  is  no  symmetry  in  y(Cos  G)  and 
therefore  ;(r ,  ',)  will  not  be  equal  to  Kr  ,-fj) .  The  symmetry 
condition,  Kqs  (3a)  and  (3b),  implies  that 

.s(r,y,x)  =  even  function  in  y  (3c) 

At  t  h r'  air  -  round  interface  :(r,T)  is  continuous  except 
at  :  -  0  ( !'e  f  6  :  !  6  0 )  ,  i  .  e  . 

(r,  )  =  tir,..)  at  z  ~  0  and  p  f  0  (4) 

air  "ri'iini 

it v  ■  r  tii'-  derivative:.  '  (  r  ,  )  are  net  continuous. 


The  problem  geometry  and  coordinate  systems  as  shown 
in  Figure  2  implies  that  when  p  is  equal  to  zero  the  angle 
X  must  also  assume  a  value  of  zero.  Therefore,  along  the 
z  axis  (p=0)  the  angle  v  does  not  vary  between  0  and  2 rr . 
This  means  that  there  is  no  x  variation  in  4>(r,f2)  at 
p=0  and  that 

=  0  for  3=0  (5) 

o  X 


\ 
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Ill  The  Finite  Element  Method 


The  finite  element  method  is  a  mathematical  and 
numerical  technique  for  approximating  the  solution  to 
a  large  class  of  problems.  Initially  ir  was  developed 
and  used  to  solve  problems  in  stress  analysis  (Ref  7:9). 
Later,  as  the  mathematical  foundation  of  the  method  was 
established  it  gained  widespread  acceptance  and  use  in 
solving  a  larger  class  of  problems. 

Finite  element?  are  an  extension  of  the  Rayleigh- 
Ritz  technique  of  first  recasting  the  problem  in  an  equi' 
valent  variational  form  and  then  seeking  a  solution  on 
the  basis  of  an  energy  minimization  principle  (Ref  8:1). 
In  the  Rayleigh-Ritz  method  a  solution  in  the  form  of 
a  linearly  independent  set  of  trial  functions  is  assumed 
These  trial  functions  must  satisfy  the  essential 
boundary  conditions.  The  approximate  or  "best"  solution 
to  the  problem  is  the  linear  combination  of  these  trial 
functions  which  maximizes  (or  minimizes)  the  variational 
principle  (functional).  If  this  linear  combination  of 
functions  is  not  an  extremum  (maximum  or  minimum)  of  the 
functional,  then  the  class  (or  space)  of  trial  functions 
is  expanded  by  the  addition  of  more  functions.  This  ex¬ 
pansion  of  the  trial  function  space  is  continued  until  a 
linear  combination  of  functions  which  is  an  extremum  of 
the  functional  is  obtained. 

'i 

The  finite  element  solution  technique  is  similar  to 
that  of  Rayleigh-Ritz.  The  only  difference  lies  in  the 
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choice  of  trial  functions.  In  finite  elements  the  problem 
domain  is  divided  into  smaller  regions  (grids)  or  elements. 
Each  trial  function  is  usually  associated  with  only  a  few 
elements.  Unlike  Rayleigh-Ritz ,  finite  elements  uses 
trial  functions  which  are  zero  over  parts  of  the  solution 
domain  (a  local  basis).  Also,  the  trial  space  (number  of 
trial  functions)  is  expanded  by  using  more  elements  (mesh 
points)  and  not  by  the  addition  of  a  new  class  of  functions. 
Because  of  these  differences  the  finite  element  method  is 
more  adaptable  towards  a  numerical  (computer)  solution  than 
Rayleigh-Ritz . 

There  are  two  basic  approaches  to  the  finite  element 
formulation  of  a  problem.  One  approach  is  to  find  the 
extremum  of  the  functional  which  originates  from  a  varia¬ 
tional  principle  and  the  calculus  of  variations  (Rayleigh- 
Ritz  with  a  local  basis).  The  other  is  by  the  method  of 
weighted  residuals.  The  method  of  x^eighted  residuals 
does  not  include  the  use  of  a  variational  principle  or 
the- calculus  of  variations.  In  some  problems  a  variational 
principle  has  not  been  developed  or  may  not  exist,  and 
therefore,  the  Rayleigh-Ritz  approach  cannot  be  used. 
However,  in  such  cases,  the  method  of  weighted  residuals 
can  be  used  to  solve  the  problem.  Therefore,  the  method 
.-,r  weighted  residuals  can  be  extended  to  a  wider  class  of 
-  ,'b loins  tn.an  Ravin ioh-R  itz  . 

fun  method  ■:  weighted  residuals  is  another  approach 
for  developing  a  sot  of  (algebraic)  problem  equations  to 
which  the  finite  element  method  can  be  applied.  There  are 


three  basic  mathematical  "recipes"  through  which  the  method 
of  weighted  residuals  can  be  developed.  These  are  the  methods 
of  least  squares,  collocation  and  Galerkin.  In  some  problems 
where  a  variational  principle  (functional)  exists  it  can  be 
shown  (Ref  9:735)  that  the  Galerkin  method  of  weighted 
residuals  is  equivalent  to  Rayleigh-Ri tz .  An  identical  set 
of  matrix  equations  and  therefore  the  same  solution  is 
achieved  by  either  method. 

In  the  following  paragraphs  a  variational  principle  for 
the  air-over-ground  problem  and  the  even  parity  form  of  the 
Boltzmann  equation  is  examined.  A  weak  form  of  the  variational 
principle  and  the  boundary  conditions  are  incorporated  within 
this  development.  Finally,  the  Galerkin  method  of  weighted 
residuals  is  discussed  and  an  equivalence  to  the  variational 
approach  for  the  air-over-ground  problem  is  established. 


Even  and  Odd  Parity  Second  Order  Forms 

In  order  that  a  variational  principle  may  be  used  the 
even  and  odd  parity  forms  of  the  anisotropic  Boltzmann 
equation  and  associated  boundary  conditions  will  be 
developed.  The  starting  point  01  this  development  is 
Eqs  (1)  and  (2).  Following  the  derivation  of  Kaplan  and 
Davis  (Ref  10:166)  and  that  oi  Wheaton  (Ref  5)  the  mono- 
enorgeiie  steady  state  transport  equation  can  be  written  in 
terms  of  ctie  -s.  vecto-  by  changing  a  to  -d  in  Eq  (1). 
this  gives 


•  V  :  (  r  .  -  ■)  ^  :  i  (  r )  ;  (  r  ,  —  ) 


f  ■ 


S(r,-M'):(r.. ')  do 


fo) 


* 


The  even  and  odd  parity  terms  will  now  be  defined  as 


Y(r . 

,Q)  = 

:  %{<K 

+ 

Kr, 

-ft)} 

(7) 

X(t; 

:  %{cp(r,S7) 

- 

d(r , 

-ft)} 

(8) 

S®  ( r . 

,«)  = 

:  %{S(r,fl) 

+ 

S(r, 

-ft)} 

(9) 

sV: 

,4)  = 

=  %{S(r,fl) 

- 

S(r, 

-ft)  1 

(10) 

C?  .  A  A 

o°(r ,a- 

■(f)  = 

C  A  A 

:  A  (  0  ( r  ,  • 

')  + 

°S(r, 

A  A 

-ft*  ft 

0) 

(ID 

au(r ,fi- 

■  r/>  = 

=  h  (  a  s  ( r  ,  ft  < 

■ft 

')  ~ 

os(?, 

A  A 

-ft*  ft 

0) 

(12) 

where 


Z'*  A 

■?(r,ro 

X(r,S) 

Sg(M) 

su(r,n) 

O  ^  /N  yN 

o°(r,.M,  ) 

U~  g 

a  (r  ) 


even  parity  fluence 
odd  parity  fluence 
even  parity  source 
odd  parity  source 

even  parity  scattering  cross-section 
odd  parity  scattering  cross-section 


Adding  Eqs  (1)  and  (6),  then  dividing  throughout  by  two 
and  using  the  above  definitions  gives 


fi*7x(r,S)  +  c  c  (  r  )  r  (  r  ,?.) 


rO  A  A  A  A  . 

o°(r.t‘li  )  e  (  r  ,  si  )  d  f. 

<  T  A. 

+  S°(  r  , ... ) 


(13) 


Using  the  derivation  of  Wheaton  (Ref  5:8),  which  is  also 
reproduced  in  Appendix  A,  the  scattering  kernel  term  in  (13) 

can  '•••  vr  1 1 r  •»» 


(r. 


)  (.  r  . 


’)  d  . 


j. 


(  r 


’)  (r  . 


’)  d. 


(14) 
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where  the  even  properties  of  the  even  parity  scattering 
cross-section  and  the  even  parity  fluence  have  been  used 
in  the  derivation  of  Eq  (14).  F.q  (13)  now  becomes 


o-7x(r,:3)  +  o  t  (  ? )  i'  (  r  ,  J  o~  (  r  ,  Q  •  S5  ')  *  ( f  ,  3  ')  df 


S"(? .3) 


(15) 


Similarly,  by  substracting  F.qs  (1)  and  (6)  and  rearranging 
the  scattering  kernel  (See  Appendix  A)  gives 


' ; ( ? ,  )  +  ’t(f)y(r,n)  = 


(  f  ,  f:  •  r;  ')  V  (  r  ,  ' )  d  ' 


n  U  ,  ^  V 

S  (r,  .) 


(16) 


Eq s  (15)  and  (16)  are  referred  to  by  Kaplan  and  Davis  (Ref  10) 
as  canonical  forms.  The  natural  boundary  condition  Eq  (2) 

Ccn  also  be  rewritten  as 


'1(?s 

C) 

V 

• 

O 

II 

/ — \ 

0 

(17) 
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e  Fig.  3).  Therefore  it 

follows  that 
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and 


xCr,^)  -  X(r,ft2) 


(21) 


Also  x(r>^)  and  T(r,P)  are  continuous  at  the  air  ground  inter¬ 
face  (when  u  ~  0)  but  their  derivatives  are  discontinuous. 

A  further  simplification  is  now  introduced  into  this 
development  by  defining  the  even  and  odd  operators, 
and  Gu  as 


G° ( r ) f ( r , P ) 
GU(r)f(r,i!) 


=  ot(r)f(r,fi)  -  f^os(r,Q'P.')£(r,n')dn' 
=  °t(r)f(r,P.)  -  ( r  0  f  (  r  ,P.  Odf. ' 


(22) 

(23) 


where  the  scattering  cross  sections  can  be  expanded  in 
spherical  (surface)  harmonics  (see  Appendix  A)  to  give 


G°(r)f(r,?n  -  ct(Of(9,P) 

-  Zv.Vf(^^m(^/nim^0f(r,^)c!.Q'  (24) 

where  the  even  parity  Legendre  expansion  cross-section  cf 
is  defined  as 

of (9)  =  for  ^  cvon 

1^0  for  li  odd  (25) 

and 


c 

c^(r)  =  Legendre  expansion  scattering  cross-sec t ion 
coeiticients  (Ref  5:2 'J) 


The  derivation 


(■■')  can  *><>  :  oun-.i  in  Append  Lx  A. 

I'Vfn 


It, 


d 


J 


The  odd  parity  scattering  cross  section  can  also  be  expanded 
to  give 


Gu(r)f(?,£)  =  at(r)£(rth 

-  iiV“G)Ws)LYL(^£(bS-)o.;;-  (26) 


where  only  odd  expansions  in  £ 


are  used  or  is  defined  as 


°oV) 


c^(r)  for  £  odd 
0  for  £  even 


(27) 


Cr  ...... 

The  G°  and  G  operators  are  self  adjoint  positive  definite 
(Ref  10:174).  Inserting  these  operators  into  Eqs  (15)  and 
(16)  they  become 


G°( r )T( r , Q)  =  S°(r,£)  -  i>VX(r,r.) 


and 


Gu(r)x(r,fi)  =  S u ( r , 3- )  -  H-VV(r ,3) 


(28) 


(29) 


Solving  for  ^(rjSl)  and  x(r ,2)  from  Eqs  (28)  and  (29)  pro¬ 
duces 


Gs(r) 


)j“]|sS(r,3)  -  3-vx(r,ro} 

'»•'-)  =  U»u(r )j  1  |su(r  ,4)  -  iJ*V’i,(r,i2)J> 


T(r,m  = 

X(r 

where  using  the  notation  of  Kaplan  and  Davis 


U  ^  11  ^  —  1 

K  (r)  =  (G  (r)}  =  inverse  of  the  operator 


G  (  r )  -  i  G  (  r  )  | 


-1 


0u ,  ~  N 
G  (r) 


inverse  of  the  operator 
CS(r) 


(30) 

(31) 


(32) 


(33) 
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A  detailed  mathematical  derivation  of  these  inverse  operators 
is  presented  in  Appendix  A.  They  are  defined  as  (Ref  11:481) 


<>  A  «*-  I  /N  1 

K°  (  r )  =  ot  ( r )  n 

♦  Vm ' °  t <  A  /  <ot  C  O <  O >  >  ')d£ ' 

(y  A  , 

where  cg(r)  is  defined  by  Eq  (25). 

1 


(34) 


KU(?)  =  (r ) 


+  E,ya^(r)/(ot(9r<(r))}\m(^/;m(^)dn  ' 


(35) 


Both  KU  and  K8  must  be  self  adjoint  positive  definite  since 
they  are  the  inverses  of  Gu  and  G8  which  are  both  positive 
definite  and  self  adjoint. 

The  functions  ^  (f2)  and  Y^  (fc)  are  the  well  known  nor¬ 
malized  spherical  (surface)  harmonics  (Ref  6:609)  which  are 
defined  as 


Vs>  '  =/2Ali  • 


(36) 


4 it  ( S.+m)  ! 


and  (ft)  being  the  complex  conjugate  of  Xm(^)  is  defined 
as 


2<’,+  1  (  0-m)  l 


4 rr  *  (£+m)  !  P£(y)e 
where  y  =  Cost)  and  P™(u)  is  the  associated  Legendre  func- 

: ions . 

Eqs  ( 3(J  and  (31)  can  now  be  written  as 


(3  7) 


■i  (  r  ,  f>  )  =  K8(r)  tss(r  ,P.)  -  il*  Vy(r  ,il)  } 


(38) 
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and 


x(r,n) 


Ku(r) {Su(r ,P) 


a-vv(r,a)} 


(39) 


Substituting  Eq  (39)  into  Eq  (28)  produces  the  second  order 
even  parity  form  of  the  Boltzmann  equation 


-S-VKu(r)a-V'F(r,«)  +  G§(r)'f(r  ,8)  =  Ss(r,ft) 

-  P* VKu(r)Su(r ,P)  (40) 


and  inserting  Eq  (39)  into  Eqs  (17)  and  (18)  gives  the  appro¬ 
priate  surface  boundary  conditions  for  Eq  (40)  as 

V(rs,a)  +  Ku(rs){S  (rs,fl)  -  a*  74f(rs  ,fl)} 

=  0  for  P*n  <  0  (41) 

and 

nr,  -  Ku(r  J{SU(r  d)  -  P.-V'^r  ,a)> 

=  0  for  >  0  (42) 


Eqs  (38)  and  (39)  give  the  second  order  odd  parity  form  of 
the  Boltzmann  equation  which  is 


a  ry  />.  , s  /s/s  j«/s  /\  /s  1 1  a  /s 

-P-VKfe(r)p.vX(r,n)  +  Gu(r)X(r  ,P.)  =  Su(r,a) 


?K#’(r)S*(r,Q) 


(43) 


Inserting  Eq  (38)  into  Eqs  (17)  and  (18)  the  surface  boundary 
condition?  for  Eq  (43)  are 


+  Ks(r){SS(r,,£)  -  d-vx(f  ,u)} 

=  0  lor  : .  • n  x  0 

(44) 

(t..,:) 

■3 

-  K*(r)  tS-:(rs,;>)  ~  r>vx(rs 

=  0  for  fi«n  >  0 

(45 ) 
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The  Variational  Principle 


Variational  principles  for  the  monoenerge t ic  transport 
equation  have  been  found.  These  principles  are  related  by 
a  series  of  transformations  (Ref  10:166).  The  functional 
whose  Euler  equations  are  the  even  and  odd  parity  second  order 
Boltzmann  equations  will  be  used  in  this  work.  Primarily 
the  even  parity  component  will  be  used  because  it  is  always 
positive,  self-adjoint  and  can  be  integrated  to  give  the 
scalar  flux  or  fluence  (Ref  12:148).  The  odd-parity  flux  which 
can  be  negative  integrates  (over  all  directions)  into  the  net 
current  (Ref  13:12).  Another  "nice"  feature  of  this  even  and 
odd  parity  formulation  is  that  it  produces  a  solution  matrix 
which  is  positive  definite  and  symmetric. 

Defining  the  inner  product  of  two  functions  as 


<f  »g> 


f?S)g(:2)dfi 


(46) 


where  *  means  the  complex  conjugate,  the  functional  which 
corresponds  to  Eqs  (40)  (41)  and  (42)  is  given  as  (Ref  10:169) 


F(u) 


yR{<*>V  u,KU(:>Vu)> 

-2<  i,S~>}dr  - 


+  <u,Cgu>  -  2<n*Vu,KuSU> 
!  6*  n  |  u“d.T)  ds 


(47) 


where  ys  represents  a  surface  integral. 

It  can  be  shown  Ref  (10:169)  that  the  Euler  equation 
(stationary  point)  of  this  functional  is  indeed  the  even  par-' tv 
second  order  term  o  f  the  Boltzmann  equation  and  that  Eqs  (41.) 
and  (42)  are  the  natural  boundary  conditions.  A  similar 
functional  for  the  odd  parity  equation  has  also  been  found. 
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In  Appendix  B  the  stationary  point  of  Eq  (4?)  is  shown  to 
be  a  minimum  and  the  weak  or  Galerkin  form  is  obtained.  This 
weak  form  is 


Li  (ft*Vn,Ku(ft  •V'T.>  +  <n 

Ocfi-Vr,  ,KUSU>  +  <n,Sg>}dr 


,Gg,F>}dr  +  ^sjhTi  I  ^*n  I  hTdftds 


h 


(48) 


where 


n  =  test  or  weight  function 

and 


T  =  trial  function 


The  natural  boundary  conditions  Eqs  (41)  and  (42)  are  incorporated 
into  the  weak  form  of  equation  (48).  However,  the  boundary 
condition  at  the  ground  interface  is  an  essential  condition. 

A  solution  to  the  air-over-ground  problem  can  be  obtained 

from  a  solution  to  Eq  (48)  and  this  essential  boundary  condition. 

The  Galerkin  or  weak  form  of  Eq  (48)  produces  a  solu¬ 
tion  matrix  which  is  positive  definite  and  symmetric.  It 
is  Symmetric  because  the  test  and  trial  functions  are  the 
same  in  the  Galerkin  solution  method.  The  matrix  is  positive 
definite  because  Ku  and  C>  are  positive  definite  and  it  is 
obvious  that  for  the  Galerkin  solution  the  term  <il  •  7n  ,K  ( f.  •  vv ) " 
is  also  positive  definite. 


The  Method  of  Weighted  Residuals 

The  methou  of  weighted  residuals  is  a  straightforward 
and  simple  prescription  for  solving  a  wide  class  of  problems. 


Unlike  Rayleigh-Ritz  it  does  not  depend  on  a  variational  prin 
ciple  or  the  calculus  of  variations.  However,  in  solving 
certain  types  of  problems  the  use  of  a  variational  principle 
or  the  Galerkin  method  of  weighted  residuals  is  equivalent 
and  they  produce  solutions  which  are  identical.  For  the  air- 
over-ground  problem  and  the  second  order  even  parity  form 
of  the  Boltzmann  equation  it  will  be  shown  that  the  Galerkin 
method  and  the  weak  form,  which  is  given  by  Eq  (43), 
are  equivalent  formulations  of  the  same  problem. 

In  the  method  of  weighted  residuals  an  approximate 
solution  which  is  a  linear  combination  of  trial  functions 
is  assumed  to  exist.  These  trial  functions  are  required  to 
satisfy  the  necessary  boundary  and  continuity  conditions. 

The  approximate  solution,  when  inserted  into  the  problem 
equation,  is  then  required  to  be  an  exact  solution  of 
the  problem  with  respect  to  several  weight  functions  (Ref 
7:106).  The  choice  of  weight  functions  determines  whether 
the  method  of  weighted  residuals  is  one  of  collocation, 
least  squares  or  Galerkin.  In  the  Galerkin  method  the 
weight  functions  are  chosen  to  be  the  same  as  the  trial 
functions . 

Applying  the  Galerkin  method  of  weighted  residuals  to 
the  second  order  form  of  the  Boltzmann  equation  is  equivalent 
to  us  inn  a  variational  principle.  In  Appendix  C  the  weak 
form,  Eq  (43),  is  derived  from  a  Galerkin  formulation  of 
this  problem.  The  natural  boundary  condition  is  incorporated 
into  this  development  and  an  equation  identical  to  Eq  (48) 


is  produced.  Therefore,  solutions  to  the  second  order  form  of 
the  Boltzmann  equation,  by  using  a  variational  principle  or 
the  Galerkin  method  of  weighted  residuals  are  equivalent. 

This  equivalency  exists  because  the  second  order  even  parity 
operator  of  Eq  (40)  is  positive  definite  and  self-adjoint. 


TV  Space-Angle  Synthesis  of  the 
Even-Parity  Anisotropic  Boltzmann  Equation 

A  space  angle  synthesis  solution  approach  has  already 
been  applied  to  the  air-over-ground  problem.  Roberds  and 
Bridgman  used  "specially  tailored"  angular  trial  functions 
to  solve  the  two  dimensional  steady  state  anisotropic 
Boltzmann  equation  Ref  (4:332).  Space  angle  synthesis  was 
applied  directly  to  Eq  (1).  and  not  to  the  second  order  form 
of  the  Boltzmann  equation.  Miller  e_t  al .  (Ref  13:12)  have 
applied  phase-space  finite  elements  directly  to  the  isotropic 
second  order  Boltzmann  equation  in  x  -  y  geometry.  Wheaton 
(Ref  5)  has  applied  phase-space  finite  elements  to  the  air- 
over-ground  problem.  However  a  space  angle  synthesis  finite 
element  approach  using  a  flux  expansion  in  bicubic  splines 
and  spherical  harmonics  has  not  been  done. 

Because  of  the  complexity  of  the  air-over-ground  problem 
a  space-angle  synthesis  finite  element  approach  seems  to  be 
justified.  This  solution  technique  might  have  several 
advantages,  some  of  which  are: 

1.  The  elimination  of  ray  effects; 

2.  The  numerical  advantages  of  finite  elements 
in  combination  with  a  space-angle  synthesis 
approach,  may  be  able  to  better  handle  the 
four-dimensionality  of  the  problem; 

3.  Anisotropic  scattering  can  be  easily  handled 
by  a  "wise  and  convenient"  choice  of  angular 
trial  f unc t ions  ; 
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4.  The  computational  effort  might  be  reduced, 
without  a  loss  in  accuracy.  It  is  expected 
that  bicubic  splines  will  not  require  a 
fine  spatial  problem  grid  (mesh  size);  and 

5.  The  boundary  conditions  and  a  first  scatter 
source  formulation  can  be  easily  handled. 

The  finite  element  space-angle  synthesis  technique  is 
merely  a  spatial  and  angular  expansion  of  the  even  parity  flux 
by  a  tensor  product  of  polynomial  functions  and  spherical 
harmonics.  In  this  work  a  tensor  product  of  bicubic  poly¬ 
nomial  splines  is  used.  This  expansion  is  the  trial  solu¬ 
tion  which  will  be  used  in  the  finite  element  method.  The 
piecewise  bicubic  spline  expansion  becomes  a  local  basis  in 
the  spatial  (p,z)  variable.  However,  the  spherical  (surface) 
harmonics  which  are  defined  throughout  the  angular  problem 
domain  form  a  global  basis.  Therefore,  this  trial  function 
expansion  has  a  dual  basis  —  a  local  basis  in  space  and  a 
global  basis  in  angle.  This  is  the  finite  element  space 
angle  synthesis  method. 

I'he  starting  point  of  this  development  is  the  weak 
form  of  the  second  order  even  parity  Boltzmann  equation 
and  essential  boundary  conditions.  The  air-over-ground  pro¬ 
blem  is  described  by  the  weak  form  of  II q  (43)  and  the 
symmetric  condition  i‘q  (20).  An  essential  boundarv  condition, 
at  the  air  ground  interface,  is  that  'l'(r,iJ)  is  continuous 
at  z  ~  0  and  f  0  (see  Fig.  _).  However,  the  derivatives 
of  7(r,k)  are  discontinuous. 


In  the  remainder  of  this  section  the  spatial  and 
angular  trial  functions  will  be  examined.  A  system  of 
coupled  equations  for  the  numerical  solution  of  the  air- 
over-ground  problem  will  be  developed.  A  first  scatter 
(collision)  source  will  be  used  in  this  development. 

The  Trial  Functions 

Because  this  is  a  four  dimensional  problem  with 
anisotropic  scattering  a  numerical  solution  technique  is 
required.  The  finite  element  formulation  of  this  problem 
lends  itself  directly  to  such  a  solution  approach.  However 
an  application  of  four  dimensional  phase-space  finite  elements 
to  Eq  (48)  will  be  very  costly  (computer  costs)  and 
inefficient  (Ref  5:33).  This  is  due  to  the  added  complexity 
of  anisotropic  scattering.  Anisotropic  scattering  increases 
the  bookkeeping  and  computational  difficulties.  A  local 
elemental  basis  in  angle  requires  that  the  scattering  contribu¬ 
tion  to  each  element  must  be  computed  on  an  element  by 
element  basis  for  all  space  and  angle  elements  within  the 
problem  domain.  Therefore  a  four-dimensional  phase-space 
finite  element  formulation  of  this  problem  is  not  a  very 
attractive  or  realistic  approach. 

A  close  examination  of  Eq  (48)  shows  that  the  problem 
operator  is  self-adjoint  positive  definite  symmetric.  This 
allow. •  _he  use  ol  standard  matrix  iterative  solution  tech¬ 
niques  such  as  Gauss-Se idol  ,  Jacobi  or  Successive  over¬ 
relaxation  (Ref  14:183).  Therefore,  a  numerical  method 
that  includes  a  finite  element  solution  might  be  feasible 


if  the  anisotropic  scattering  contributions  can  be  effectively 
dealt  with.  A  space-angle  synthesis  finite  element  development 
with  a  local  spline  basis  in  space  and  a  global  spherical 
harmonic  basis  in  angle  appears  to  meet  this  requirement. 

A  phase-space  finite  element  problem  formulation  which 
eliminated  the  characteristic  lines  of  the  hyperbolic  discrete 
ordinate  equations  has  been  effective  in  mitigating  ray 
effects  ( Re i  3:205).  The  well  known  and  double-P^,  equations 
of  nuclear  reactor  physics  have  inherent  elliptic  features 
which  eliminate  ray  effects.  Therefore  space-angle  synthesis 
using  spherical  harmonics  seems  to  represent  an  approach 
which  will  eliminate  ray  effects. 


The  trial  function  expansion  which  will  be  used  in  this 
development  is 


T.Z 


1R  +L  +  f 


i(  s ,  P 


>  P  >  X )  -  V  Y  y*  y  Aiz,irBiz^Z')'Bir^°'>»mi^,x^ 

L — i  L — i  Z _ <  L _ ,  t,m 


iz=l  ir=l  i=0  m=— i 


(4b) 


where  c,z  is  the  spacial  coordinate  dependence  in  cylindrical 
geometry  and  u,X  represents  the  2  angular  variable  in  Fig.  2. 


V  (  z  >  P  »  V  ,  X ) 


A- 

iz  ,  lr 

^  jHl 


BLz(z) 


3lr(0 


even  parity  fluence  at  position  r,z  and  in 
direction  y,x> 

flux  expansion  or  mixing  coefficients, 

cubic  polynomial  spline  in  the  in¬ 
variable  (z-spline), 

cubic  polynomial  u-spline, 

spherical  harmonic  function,  Eq  (51). 
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Figure  4.  Cubic  Spline  With 

Evenly  Spaced  Nodes 


iz,  ir,  £  and  m  are  the  trial  function  summation  indices  and 

IZ  =  total  number  of  z-splines, 

IR  =  total  number  of  p-splines, 

L  =  degree  of  the  spherical  harmonic  expansion. 

The  definition  of  a  cubic  polynomial  splines  (Ref  15:89) 
with  evenly  spaced  nodes  (knots)  is 


^i+ 


xi+  <X«Xi^ 


(xi+  -x)3-4(xi  +  i-x)\ 


x-«:x4x^+ 


j  c  >:  •  .  )’-V  X  •  +  _  -x  )'46(  X  X  )  ,  x  ;  _  , -SX*  X  • 

1  Oc«x-_ 

r:r.m  oi  t  a  l  s  -  r  i.  :  fune  t  i  on  is  shown  i  n  Fie.  4.  l"n<? 


Y£m(y’x)  =  C£mP£m(u)elmX  =  SnP£m(lJ)  {cos(mx>  +  lsin<mx>} 


(51) 


where 


C£m 


(£  ~  K  ! 

(£  +  m) ! 


(52) 


The  trial  function  expansion,  Eq  (49)  can  be  made  to  satisfy 
the  symmetry  of  Eq  (20).  This  is  a  symmetric  condition  in  x 
which  directly  implies  that  the  solution  must  also  be  an  even 
function  in  the  variable  X.  Therefore,  the  angular  trial 
function  expansion  of  Eq  (49)  must  also  be  even  in  y •  Dropping 
the  isin(m\)  term  from  Eq  (51)  and  substituting  into  Eq  (49) 
gives 


'•'(z.P^.X) 


LR 

V 

/ 


iz=l  ir=l  l-'O  m=0 


V"  V  A-  •  (z)' 

f  J  vz,  lr  rzv  ' 

£  ,m 


B-  (P)*Q, 
irv  '  ^xvi 


(33) 


where 

C4m  =  C;mP.'.m(u)cos(mx)  (54) 

ana  by  using  the  or thonorma 1  properties  ol  spherical  harmonics 
the  m  index  begins  at  zero  instead  of  -£  (see  Appendix  D) . 

The  essential  boundary  condition  at  the  air  ground 
interface  must  also  be  applied  to  Eq  (53).  The  fluence 
coat  inu :  ty  requirements  can  be  satis  t  ioa  by  this*  expansion  in 
bicunic  polynomial  splines  and  spherical  harmonics.  Both 
•  ; f  tiies  functions  ar  ■  continuous  throughout  the  pro D lorn 


dorna  i  n  . 


The  cubic  splines  are  also  twice  continuously 


differentiable  but  the  spatial  z-derivative  of  the  solution 
fluence  is  not  continuous  at  the  ground  interface.  However, 
the  z-splines  can  be  modified  to  have  discontinuous  deri¬ 
vatives  at  this  interface  (Ref  16).  A  Double-P.,  or  Yvon's 
method  (Ref  6:163)  can  also  be  used  to  accomodate  the  fluence 
discontinuity  at  p  =  0.  In  this  development  the  air-ground 
interface  will  not  be  included  in  the  problem  domain,  and 
therefore,  this  interface  boundary  condition  will  not  be 
enforced . 

Since  a  Galerkin  method  is  being  used  the  test  or  weight 
functions  are 


■l(z,p,p,x)  =  B  j  z  (  z)  Bjr(p)Qkn  (55) 

Substituting  Eqs  (55)  and  (53)  into  Eq  (48)  produces  Eq 

(58)  where  for  simplicity  the  p,z,p,x  dependencies  are  omitted 

and 


and 


i  =  Biz(z>Brr<p> 


(56) 


3  •  =  ,3.  ( a  1 3  •  <  ) 

■  1 1  •/  \  —  /  .  y.  \  •  / 


(57) 


i  K.n  ‘ 


i  'JCm/ 


...d  ds! 


(B-t)  '  ,RU5U 


R  -  ’  '•‘-'j’'  -o  ^  J 


(5°) 
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This  is  a  system  of  coupled  algebraic  equations  where  the 


unknown  quantities  are  the 
can  also  be  written  as 


Aijrn  mixin§ 


coefficients . 


Eq  (58) 


NxN  Nxl  Nxl 


where 


(59) 


N  =  IZ • IR 4 (L+l) • (L+2)/2  (60) 

K  =  Coefficient  or  stiffness  matrix,  where 
each  element  is  a  summation  of  terms  in 
the  square  brackets  of  Eq  (58) 

A  =  A.-0  mixing  coefficient  vector 

1.  )  A/  m  v 

S  =  Source  vector  which  is  the  right  hand  side 
of  Eq  (58) 


The  A  £ ; p  m  coefficients  of  Eq  (53)  will  be  obtained  from 
a  computer  solution  of  Eq  (59) «  These  mixing  coefficients 
can  then  be  substituted  into  Eq  (53)  to  give  the  even  parity 
fluence,  i (z , p , p, x) .  In  a  solution  of  Eq  (59)  the  elements  of 
the  K  matrix  and  S  vector  must  be  computed.  This  computation 
involves  an  evaluation  and  summation  of  the  individual  expanded 
terms  of  Eq  (58).  This  expansion  is  carried  out  in  Appendix 
1/ 

The  directional  gradient  operator  in  cylindrical  geometry 
is  defined  as  (Ref  6:59) 

-W;  =  V  l-T~cosv  2i£±)  -  1  3{)Y  I-r  s  i  n \  }  +  u3  :•  (61) 


i'his  is  the  conservative  form  ot  the  directional  derivative 
in  two  dimensional  (p,z)  geometry  with  azimuthal  symmetry. 
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The  Ku  and  operators  have  been  defined  in  Eqs  (24)  and 
(35)  of  Chapter  III.  The  scalar  product  of  the  velocity  and 
normal  vectors  is  given  in  Appendix  E  as 


o*n 


for  the  horizontal  outer  surfaces  (top  or 
bottom)  of  the  problem  cylinder,  and  as 

for  the  vertical  surface  (side)  of  the 
cylinder  (62a) 


where 


:  *n_ 


;•  on  the  top  surface,  and 
-p  on  the  bottom  surface 


(62b) 


and 


=Vl  ~U2  cosx  (62c) 

P 

The  normal  unit  vectors  h^  and  h^  are  shown  in  Fig.  4,  Appendix  E. 

Expanding  the  expressions  in  Eq  (58)  produces  an  integral- 
differential  equation  which  has  twenty-eight  terms  (see  Appendix 
E) .  These  terms,  except  for  the  source  terms,  can  be  easily  se¬ 
parated  into  a  product  of  z,  p,  \i  and  y  integrals.  This  is  an 
integral  separation  of  variables  which  is  a  direct  result  of  Eq . 
(53);  where,  it  is  assumed  that  the  solution  can  be  expressed  in 
a  form  where  the  dependent  variables  are  separable.  This  separa¬ 
tion  property  simplifies  the  individual  integrals  which  have  to 
be  evaluated.  It  allows  for  the  evaluation  of  only  single  inte¬ 
grals  and  not  the  more  complicated  double,  triple  or  quadruple 
integrals.  By  this  separation  of  variables  it  may  be  possible 
to  integrate  most  of  these  single  integrals  analytically  and 
thus  avoid  a  numerical  integration  process. 

The  Spherical  Harmonic  Integrals.  The  use  of  a  spherical 
harmonic  angular  trial  function  expansion  was  motivated  by  six 
basic  cons  i  ue  rn  1 1  one-  : 


a  a 


1.  Because  of  the  global  nature  of  these  functions  the 
computational  effort  will  be  substantially  reduced. 

2.  Spherical  harmonics  are  well-known  functions  with 
orthonormal  and  symmetric  (odd,  even)  properties. 

3.  The  scattering  cross  sections  are  usually  expanded  in 
spherical  harmonics  (see  Appendix  A). 

4.  Spherical  harmomics  will  produce  a  system  of  equations 
which  are  elliptic  and  invariant  under  continuous 
coordinate  rotations  (Ref  1:362). 

5.  An  analytic  or  closed  form  evaluation  of  the  angular 
integrals  might  be  possible. 

6.  The  even  parity  angular  fluence,  Eq  (53),  can  be 
easily  integrated  to  give  the  total  particle  fluence. 

This  integration  is  carried  out  in  Appendix  I. 

The  term-by-term  expansion  of  Eq  (58)  has  been  partly  carried 
out  in  Appendix  E.  The  resulting  angle  integrals  are  only  de¬ 
pendent  on  the  degree  of  the  spherical  harmonic  trial  function 
expansion  which  is  used.  They  are  not  dependent  on  the  problem 
parameters  and  therefore  they  can  be  independently  evaluated. 

They  can  be  evaluated  once,  and  thereafter,  used  as  a  part  of 
the  problem  input  data. 

Three  approaches  were  pursued  in  an  attempt  to  evaluate  the 
angle  integrals  which  are  produced  by  this  expansion.  The  first 
approach  was  to  use  the  orthogonal  properties  of  the  associated 
Legendre  functions  and  the  well-known  properties  of  sines  and 
cosines  to  analytically  evaluate  these  integrals.  However,  be¬ 
cause  of  their  complicated  nature  (see  Appendix  F)  a  closed  lorm 
analvt  i  c  in  to  r  a  t  i  on  was  not  easily  obtained  for  most  o!  them. 

i'iie  second  approach  was  to  use  a  computer  routine  that  can 

•  *va  lua  i  e  they'  integrals  in  a  symbolic  or  algebraic  sense.  Sues 

a  routine  will  transform  the  integrals  into  algebraic  expressions, 
iii-  'o-o  Ha  evs  sin  wnich  was  written  bv  the  Massachu.se  i  t  s  !  ns  ‘  :  - 

•  g*  > :  !  ■  ;  !■■■  g.iLh  Arouu  was  us-si  in  this  aLt  i.mt.  Pouvve- , 

J3 


because  of  the  time  constraint  on  this  research  project  and  the 
need  to  learn  a  new  programming  language  this  approach  was  aban¬ 
doned  . 

Finally,  it  was  decided  to  evaluate  these  integrals  by 
a  numerical  integration  technique.  Because  it  is  possible  to 
separate  the  integration  variables,  the  computational  effort 
can  be  substantially  reduced  by  using  a  single  (one  variable) 
integration  routine.  So  as  to  further  reduce  the  computational 
effort,  Eq  (58)  was  completely  expanded  and  twenty  distinct 
angle  integrals  were  identified.  These  integrals  can  be  found 
in  Appendix  F.  A  ten  point  Newton-Cotes  single  integration 
routine  was  used  to  evaluate  them.  They  were  evaluated  for 
each  combination  of  the  m,  ,  k  and  n  trial  and  weight  function 
subscripts. 

Bicubic  Polynoni ai  Splines.  The  spatial  dependence 
of  the  particle  fluence  in  the  air-over-ground  problem  is 
approximated  by  a  product  of  cubic  splines.  The  use  of 
cubic  polynomial  splines  in  the  trial  function  expansion  of 
Eq  (53)  requires  the  formation  of  a  tensor  product  space, 
ibis  sou cc  is  r.iauo  u;  of  bicubic  polynomial  splines  which 
are  products  oi  o  u;ui  z-splines  on  a  rectangular  grid 
(Kef  15:131).  The  exact  shape  of  these  bicubic  splines  are 
obtained  iron  a  variational  principle  or  the  equivalent  method 
of  i  .no  ■  r-’si  iiia  1  •»  (!.<•.,  a  solution  of  Eq  (58)>l. 

jioiiiiv  p<  ■  1  vnoiii .  a  1  splines  uro  being  used  for  the 

1.  i'li.  ■.  :ir  ■  e.  uu wise  continuous  and  form  a 
>  '  -on  t  :  he  in  t  ra  I 

■  ■  ■  :  ‘  “  1  i  '  ->  . 
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This  reduces  the  number  ol  integrals  which  must  be 
evaluated  and  also  produces  a  sparse  and  banded 
coefficient  matrix. 

2.  A  separation  of  the  p  and  z  integration  variables 
is  possible. 

3.  Third  degree  polynomials  such  as  cubic  splines  have 
a  faster  rate  o':  convergence  than  those  of  lower 
degree.  Cubic  splines  are  also  twice  continuously 
differentiable  and  thus  they  are  very  smooth  func-N. 
tions.  For  a  given  problem  mesh  size  splines  will 
produce  a  coefficient  (stiffness)  matrix  that  is 
smaller  but  less  sparse  than  hermites  or  lagrange 
polynomials . 

The  expansion  of  Eq  (58)  with  a  trial  function  of  bicubic 
splines  and  spherical  harmonics  is  carried  out  in  Appendix  E. 

A  further  expansion  of  these  equations  and  a  separation  of  the 
variables  of  integration  produced  seventeen  distinct  c  and  z 
integrals.  These  integrals  which  include  the  space  source 
integrals  are  listed  in  Appendix  G.  The  source  integrals  are 
derived  from  an  interpolation  of  the  first  scatter  source  over 
the  entire  spatial  problem  domain. 

The  Source  Trms.  A  numerical  solution  to  the  air-over¬ 
ground  problem  and  the  second  order  Boltzmann  equation  requires 
that  the  source  terms  (right  hand  side  of  Eq  (58))  must  be 
evaluated.  Those  terms  form  the  individual  elements  of  the 
problem  source  vector  ir,  Eq  (59).  The  even  and  odd  parity 
sources  S"  amd  S  will  be  defined  as  the  first  scatter  or 
collision  even  and  odd  parity  sources.  The  first  scatter 
source  S(r,.:)  is  the  number  density  of  particles  which  leave 
the  burst  point  and  undergo  only  one  collision  before  being 

.  .  ,  '  .  ,  A  . 

scattered  into  direction  '  at  position  r.  Streaming  neutrons 
which  leave  the  burst  point  and  do  not  collide  before  reaching 
position  (r,  .)  are  not  included  in  the  collision  source. 
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The  use  of  a  first  scatter  source  makes  the  air-over¬ 
ground  problem  more  isotropic.  It  removes  the  strongly 
anisotropic  streaming  particles  from  being  a  part  of  the 
problem  source.  Therefore,  the  solution  fluence  of  Eq  (58) 
will  be  the  scattered  even  parity  fluence  'i  s(r  ,Q)  and  not 
the  total  even  parity  fluence  V  (r,fi).  The  total  even  parity 
fluence  can  be  defined  as 

Vt(r,fl)  =  ¥s(r,fl)  +  ^(r.fl)  (63) 

where 

V  ,  ( r  , A )  =  streaming  uncollided  particles  at 
position  ( r , A)  . 

A  precise  mathematical  definition  of  the  S®  and  Su 
sources  will  now  be  developed.  Also  a  source  interpolation 
procedure  will  be  outlined.  This  source  interpolation  is 
used  in  order  to  simplify  the  source  integrals  of  Appendix 
E  ( E- 40  to  E-46) • 

The  First  Scatter  Sour  ^ .  The  even  and  odd  parity 
sources  have  been  defined  as 

SU(r,~)  =  1;  { S  (  r ,  A )  -  S(r,-A)}  (lQ) 

SS(r,..)'=  ^{S  (r,£)  +  S  (  r  ,  -  A)  }  (9) 

JJ  O  .  A  A 

IL  S  and  S°  are  first  scatter  source  densities  then  S(r,.;) 
and  S ( r , -  . )  must  also  be  defined  as  first  scatter  source 
part  ic  ies/uni  l  volume.  If  a  position  (•',•_)  in  the  problem 
domain  is  chosen  then  a  unit  vector  from  the  burst  point 


(0,zb)  can  be  detined  as 


Figure  5.  First  Scatter  (Collision) 
Source  Direction  Vectors 


pep  +  (z-zb)ez 
(p2  +  (z-zb) 2 


(64) 


Fig.  5  shows  the  direction  vectors  of  this  first  scatter 
(collision)  source,  fi  is  the  direction  that  all  streaming 
(uncollided)  particles  have  at  point  (p,z). 

By  definition  only  particles  which  are  streaming 
radially  outward  from  the  burst  point  can  be  included  in 
the  direct  fluence.  Therefore  the  direct  fluence  at  point 
(p,z)  is  in  the  fi'  direction  and 


can  be  written  as 


(65) 


$j(p  >z  >p‘ ') 


=  _L,  exp{'fS  o  (z)dsl 

<HTT  J0  L 


where 


{p2  +  (z-zb)2P 


(66) 


and  Jd s  means  that  the  integration  is  carried  out  along  the 
path  s  (see  Fig.  5).  Also 


Y  =  Total  particle  yield  of  the  nuclear 
explosion  at  the  burst  point. 


at(z)  = 


ct(0)e“z/sh  for  z  >  0 
a  (ground)  for  z  <  0 


sh  =  atmospheric  scale  height 


(67) 


The  term  t( z)ds  is  the  average  number  of  collisions 
which  a  particle  undergoes  in  traveling  from  the  burst  point 
(o,zb)  to  point  (p,z).  From  Fig.  5  the  distance  s  can  also 
be  written  as 


s  =  (z-zb)  /’id  (68) 

and  therefore  by  changing  variables 

ds  =  dz/ud  (69) 


where  ud  is  a  function  of  o  and  z  (but  constant  along  a  path 
lengh  8) 


yd(^,z)  -  cos 


(z-zb)  /  s  =  (2-20)/  (p2  ( ::-zb):  }  (70) 


The  integral  term  of  Eq  (65)  can  now  be  written  for  z  >  0  as 

U«»fr 


(71) 
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and  finally  as 


/  \  0t^°^r  -zb/sh  ^-z/sh-i 

t(p,z)  =  — {e  ~e  }• 


sh 


=  {o^zb)  -  ot(z)}*sh 

HcT 


(72) 


From  the  above  derivation  it  follows  that 


T(p,z)  =  l 


(ot(zb)  -  at(z)}‘sh/yd  for  z  >  0 

(a  (zb)  -  at(0)}*sh/yd  -  atz/yd  for  z  <  0 


(73) 


also 


4>j(D,z,n')  = 


4  71$ 


exp  ( -t ( p  ,  z ) ) 


(74: 


Mote  that  q^(  p  ,  z  ,  S7  ")  is  only  a  function  of  p  and  z. 

The  first  scatter  source  at  (p,z)  and  with  direction 

A 

Q  are  those  particles  which  undergo  their  first  collision 

.  A  A 

at  (p,z)  and  are  scattered  from  direction  i 1'  to 
Therefore  the  first  scatter  source  can  now  be  defined  as 


S(p,z,s1)  =  as(z  ')bd(P  ,z  ,P-  ') 


(75 


S  ,  A 

where  a  is  not  a  function  of  Si'  but  of  the  scattering  angle 

o  c 

a  N 

£  *  i)  '  ■'  )  and  z.  From  Fig.  5  and  Fig.  2  '  is  defined  by 
yd  and  x  =  0  i.e.  Si'  =  (p‘',x‘*)  where  ]i  '  =  yd  and  \'  =  0. 

By  use  of  the  addition  theorem  it  is  shown  in  Appendix 
H  that  Eq  (75)  can  be  written  as 


S(- 


^  -z/shV'V' 

7  > 

L _ 


•  m 


(y)F;in(Fd)cosmX 


and  that  the  even  and  odd  parity  first  scatter  sources  are 


L  Z 


S8(P,2,S)  =<.d(P,Z,3')e-z/sh^yO=(0)C[n(l*(-l),'-m}Pta(l.d)Ptm(u)cosmx 

1=0  m*=0  (77) 


and 


L 


ii  ~  ^  — z  /  sh\  *  \ 

S  (o,z,w)  -  y-d(P,z,!i  )e  ^  J  v  ^idn^'^m' 

£=0  rff“=0 


°o (0)C“m{l-(-l)'  }Pf m(yd)P5m(p)cosrnx 

(78) 


where  m"  means  that  all  terms  with  a  m  =  0  subscript  must  be 
divided  by  two,  and  cosmx  should  be  interpreted  as  cosine  (my) 
Source  Interpolation.  Because  of  the  complicated  nature 
of  the  source  expressions  ;  Eqs  (77)  and  (7S),  and  the  need  to 
integrate  the  source  terms  of  Appendix  E,  a  spatial  source 
interpolation  will  be  used.  This  interpolation  which 
simplifies  the  source  integrals  is  necessary  if  a  very 
tedious  (double  or  quadruple)  integration  is  to  be  avoided. 

By  this  interpolation  process  the  source  terms  of  Appendix  E 
can  all  be  separated  into  a  product  of  single  integrals. 

It  is  important  to  note  that  pd  which  is  given  by  Eq  (70) 

is  a  function  of  p  and  z  and  therefore  P„  (pd)  is  also  a 

A,m 

function  of  o  and  z.  Furthermore  ip.,(p,z,  ')  of  Eq  (74)  is 

a  function  of  p  and  z.  Beginning  with  Eqs  (77)  and  (78) 
they  can  be  rewritten  as 


Sg(e 


L  p 


)P  (p)cosmx 


(79) 


S“(o,z,P.)  = 


af(0)CT,  (l-(-l) >Aflm(P,z)Pflm  (y)cosmx 


£=U  ni- ■'=(.» 


where 


Anm(p,z)  =  )P„  (yd)e 


-z/sh 


A  spatial  (p,z)  interpolation  of  the  even  and  odd 
parity  sources,  Eqs  (79)  and  (80),  is  therefore  an  inter¬ 
polation  of  (p,z).  In  this  project  these  first  scatter 
second  order  sources  were  interpolated  by  a  combination  of 
piecewise  bi~linear  Lagrange  polynomial  functions.  Speci¬ 
fically,  S®  was  approximated  by  a  tensor  product  of  linear 
Lagrange  polynomials  as  follows 


::z  nr 

\  \  <y  ' 

yP,z,  )  - ^  S3(pj,ZiP 


>)Hj(p)Hi(z) 


where 


i=l  j=l 


S®(P  •  ,Zj;  ,1;)  =  the  even  parity  source,  Eq  (79) 
evaluated  at  the  spacial  nodes 
(n ; ,Zi ) 

J 

NZ  ~  total  number  of  z-nodes 
NR  =  total  number  of  R-nodes 
l!(p)  =  Q-linenr  Lagrange  polyn<i:ial 
L(.:)  -  ^-linear  Lagrange  polynomial 
;Sj  linear  polynomials  are  delined  as 


'  < 


■'  •  ■  .  n  i  U  .  i' \  —  f  — 


■yl~l 


1 


F  i  g  u  r  e  6  . 


A  Linear  Lagrange  Polynomial  Function. 


The  product  Hj(c)-H^(z)  of  Eq  (82)  forms  a  tensor  product 
space  on  a  rectangular  grid,  in  the  p,  z  plane  (Ref  15:129) 
Substituting  Eq  (82)  for  S8  (  ,' j  ,  z  i  .  T. )  in  Eq  (77)  gives 

L 


S8( 


>  ■  >:  m 


(o)crm a 

xm 


KZ  NR 


x -0  iii'-u  + 


v,m 


- ) cosmv^  ^  2., m( P  j  > z j_ )H j  (P )Hj_ ( z )j 


i=l  j=l  (84) 

Similarly  the  odd  parity  source  can  also  be  expressed  as 


Su(.  £  [^(O)C^ 

-(-1) 


.,{1 


=U  -U 


xm 


NZ  NR 


F ,  (u)cosnv 


an 


:.=i  i=i 


(35) 


•at  -d  L n i  o  the  -ource  terms 
-  ub.;  t  :  tut:  on  and  a 


j'-aruti..  ;  t  >  :rui  •  .  ••  v.-.  r  i  .teles  produced  six  distinct 

space  inteera  1  s  .  I  iiese  space  i  n  t  ogra  l  s  are  listed  in  Appendix 


G.  The  source  angle  integrals  have  been  included  in  those 
integrals  which  are  presented  in  Appendix  F. 


V  Computer  Implementation  and  Results 


A  numerical  solution  of  the  air-over-ground  problem 
can  be  obtained  from  a  computer  solution  of  the  system  of 
coupled  algebraic  equations  which  are  represented  by  Eq 

(58) .  These  equations  can  be  written  in  the  matrix  form 

of  Eq  (59)  and  through  a  direct  inversion  or  some  other 
matrix  solution  process  the  A.  ^  mixing  coefficients 

can  be  found. 

The  computer  solution  to  this  problem  was  accomplished 
by  a  two  step  process.  Each  element  of  the  stiffness  matrix 
and  source  vector  was  computed  and  assembled  in  the  matrix 
form  of  Eq  (59).  Then  the  mixing  coefficients  were  computed 
by  the  iterative  matrix  solution  method  of  successive  over¬ 
relaxation.  An  indirect  iterative  matrix  solution  method 
is  possible  because  Eq  (58)  and  the  stiffness  matrix  of  Eq 

(59)  is  symmetric  positive  definite.  A  computer  program 
which  assembles  the  problem  matrices  and  computes  the 
mixing  coefficients  and  particle  fluences  was  written. 

This  program  which  is  written  in  FORTRAN  V  is  listed  in 
Appendix  J. 

Using  a  ten  point  Newton-Cotes  numerical  integration 
routine,  each  of  the  thirty-seven  integrals  in  Appendix 
F  and  C  were  evaluated.  The  angular  integrals  were  evaluate 
for  each  ..m,  kn  combination  oi  the  trial  ana  weight  function 
subscripts .  Selected  products  ol  these  integrals  were  then 
used  to  generate  each  of  the  twenty-eight  terms  (E-i9  to 
E-Uh)  of  Appendix  E.  Following  tin:  prescription  of  Appendix 


E  the  first  twenty-one  terms  were  then  added  (or  sub¬ 
tracted)  to  produce  the  elements  of  the  stiffness  matrix. 

The  next  seven  terms  gave  the  elements  of  the  source  vector. 

Writing  the  synthesized  Boltzmann  equation,  Eq  (58),  in 
operator  notation  as 


L(  '  t'  •  •  )  n  ■  •  =  Sn  ■ 

v  i  x  ,  i  r  '  j  z  ,  ]  r  J  z  ,  |  r 

•  ,m  k,n  k ,  a 

where  T  and  "  are  defined  by  Eqs  (53)  and  (55)  and 

L  (  i'  -  •  )  ’1  ■  .  =  Left  hand  side  of  Eq  (58) 

1  iz.ir;  jz.jr, 

£  ,  m  k ,  n 

Sn  •  =  Rivht  hand  side  of  Eq  (58) 

j  z ,  i  r  -  n 

k ,  n 

the  K  element  of  the  stiffness  matrix  is 

p,q 


(B6) 


Kp,„  ■  I-^i2(OBir(p)Q;m}-5j2(Z)Bjr(P)Qkn 


(37) 


where 


q  =  (m+1)  +  i(f+l)  +  (Lmax  +  l)(Lmax  +  2)‘{(ir-l)  +lR(iz-l)}  (88) 


and 

p  =  (a*!)  -  k(k+ 1 )  *■  iLr.iax  •«-  I)(linax  +  d)«i(jr-I)  +  JR(jz-l)}  v.  '-9) 


The  corresponding  source  vector  element  S,  is  given  bv 

Sp  -  S-Bj2(Z)Bjr(P)Qkn 

whore  p  is  given  by  Eq  (89)  and 


■*  ) 


(90) 


Lmax  =  degree  of  the  spherical  harmonic  expansion 
JR=IR  =  total  number  of  p-splines 


p  =  the  row  index  of  the  problem  (K)  matrix 
q  =  column  index  of  matrix  K 

B(z),  B(:j)  and  Q^m  are  defined  in  Chapter  IV  (Eqs  (50)  and 
(54)).  iz,  ir,  jr,  jz,  Z ,  m,  k,  and  n  are  the  trial  and 
weight  function  expansion  subscripts  where 

f.,k  =  0  to  Lmax 

m  =  Q  to  £ , 
n  =  0  to  k 
iz,jz  =  1  to  IZ 
ir,jr  =  1  to  IR 

and 

IZ  =  total  number  of  z-splines 
IR  =  total  number  of  p-splines 

I'sing  the  notation  of  Eqs  (87),  (88)  and  (S9)  the  K-problem 
matrix  and  S-source  vector  can  be  easily  assembled  in  the 
following  do  loop. 


DO 

10 

jz  =  1, 

IZ 

DO 

10 

jr  =  1 , 

IR 

DO 

10 

k  =  o, 

Lmax 

DO 

10 

n  =  0  , 

k 

fi  =  S-.v 
p 

n . 

1  z .  1 

DO 

10 

i  z  =  L  . 

lZ 

DO 

10 

ir  =  1, 

IR 

DO 

10 

:i  --  0, 

Lmax 

DO 

10 

m  =  0 , 
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p,q 


LC?, 


iz,ir' 


n  ; , 


i  J  r 
k,n 


10  continue 


The  assembled  coefficient  (K)  matrix  which  results  from 
Eq  (58)  and  a  bicubic  spline  spherical  harmonic  trial 
function  expansion  is  sparse  and  symmetric.  Because  of  this 
symmetry  and  sparseness  the  coefficient  matrix  can  be  stored 
within  the  computer  by  special  storage  schemes  (Ref  14:70). 
These  sparse  and  symmetric  schemes  will  greatly  reduce  the 
computer  storage  requirements.  As  an  example,  in  a  trial 
function  expansion  where  JR  =  IR  =  8  and  Lmax  =  2  the 
problem  coefficient  matrix  is  a  384  x  384  square  symmetric 
matrix  with  many  elements  which  are  zero.  Therefore  if  this 
entire  matrix  is  stored  within  the  computer  it  will  req.  e 
147456  separate  storage  locations.  This  much  core 
storage  is  already  beyond  the  capacity  of  most  computers.  How¬ 
ever,  by  using  a  sparse  and  symmetric  storage  scheme  this 
matrix  can  be  reduced  to  one  with  less  than  73728  elements. 

For  large  trial  function  expansions  (JR  =  IR  =  50,  Lmax  =  3), 
special  auxiliary  storage  and  solution  techniques  will  be 
necessary. 

This  entire  problem  (coefficient  and  source  matrices) 
was  assembled  on  a  CDC  6600  computer  at  the  Air  Force 
Institute  of  Technology.  The  384  x  384  problem  matrix 
is  too  large  to  be  stored  in  core  memory  unless  a  sparse 
symmetric  storage  mode  is  used.  Because  some  of  the  (p) 
integrals  in  Appendix  G  are  discontinuous  (infinite)  at 
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P  =  0  it  was  necessary  to  use  a  lower  p-integration  limit 
of  1.0  E-8.  Also,  since  the  first  scatter  sources  of  Eqs 
(77)  and  (78)  are  undefined  at  the  burst  point  none  of  the 
problem  nodes  can  be  located  there  (see  Fig.  5). 

Results 

The  computer  routines  which  are  listed  in  Appendix  J 
were  used  to  produce  a  numerical  solution  to  the  air-over¬ 
ground  problem.  These  routines  demonstrate  the  feasibility 
of  using  FESAS  to  produce  a  computer  solution  to  the  two- 
dimensional  steady  state  anisotropic  Boltzmann  equation. 

This  computer  program  has  not  been  fully  developed,  refined  or 
debugged  and  therefore  the  accuracy  of  the  results  has  not 
been  evaluated.  These  results  are  presented  in  an  attempt  to 
further  show  that  FESAS  is  a  viable  numerical  solution 
toehn  ique . 

The  problem  domain  is  a  cylinder  which  sits  on  the  sur¬ 
face  of  the  earth  (see  Fig.  5).  However,  the  air-ground 
interface  is  not  included  in  the  problem  domain  and  there¬ 
fore  all  ground  effects  are  ignored.  The  following  problem 
parameters  were  used 

V.’eapon  yield  =  l.OE+23  neutrons 

Cylinder  height  =  .4km 

Cylinder  radius  =  .4km 

Burst  height  =  .12  km 

Total  cross  section  ( g ( ( 0 ) )  -  15.0  km  ' 


Table  I 


Legendre  Expansion  Coefficients  which  were  used 
in  a  Numerical  Solution  of  the  Air-Over-Ground  Problem 


Expansion  subscript  £ 

Legendre 

expansion  coefficients 

Q 

^09 

u 

0 

i 

10.0 

.0 

! 

!  1 

0.0 

2.5 

The  c ros s-sec t ions  in  Table  I  were  arbitrarilv  chosen  and  they 
do  not  represent  the  actual  values  for  air  at  sea  level.  A 
relative  convergence  criteria  (.001)  which  is  accurate  to  three 
significant  figures  was  used. 

Ihe  program  execution  times  varied  with  the  degree  of  the 
spherical  harmonic  trial  function  expansion  and  the  problem 
mesh  (grid)  size.  The  entire  problem  matrices  were  stored 
within  core  memory  and  by  trial  and  error  it  was  determined 
that  a  relaxation  parameter  of  1.7  gave  the  fastest  convergence 
rate.  However,  as  more  trial  functions  were  used  and  the 
system  or.  equations  anu  matrices  grew  larger  the  rate  of 
convergence  decreased.  In  Table  II  the  program  execution  times 
and  the  number  of  iterations  to  convergence  are  listed  for 
^-ohl.-m  mesh  sixes.  These  execution  times  and  con¬ 
vergence  rates  can  be  subsLant  in  l  ly  reduc'd  by  rewriting  or 
"tod  !  t  v  i  ng  i  he  nre-Tai’i  reu  L  i  ne---  ;  xuoend  i  . 


FLUENCE^ (NEUTR0NS^KMhk2) 


PRRTICLE (NEUTRON)  FLUENCES. 


RFIDIUS(KM) 

Figure  7.  Computed  Fluences  as  a  Function  of  Radius  and 
showing  a  variation  with  the  Problem  Spatial 
Mesh  Size. 
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PARTICLE (NEUTRON)  FLUENCES. 


0.2  0.3 

RADIUS (KM) 


Figure  8.  Computed  Fluences  as  a  Function  of  Radius  and 

showing  a  variation  with  the  Degree  of  the  Angular 
Spherical  Harmonic  Trial  Function  Expansion. 
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In  Figures  7  and  8  fluence  values  as  a  function  of  radius  are 
presented.  These  values  are  representative  of  an  altitude  of 
200  metres  and  the  aforementioned  problem  parameters.  Figure 
7  shows  a  variation  of  the  solution  with  spatial  grid  (mesh) 
sizes  whereas  Figure  8  shows  a  variation  with  the  degree  of  the 
spherical  harmonic  trial  function  expansion.  More  detailed 
numerical  results  can  be  found  in  Appendix  K. 


VI  Conclusions  and  Recommendations 


Conclusions 

A  finite  element  space-angle  synthesis  (FESAS)  solution 
of  the  steady  state  anisotropic  Boltzmann  equation  in  two- 
dimensional  cylindrical  geometry  has  been  presented.  In  this 
presentation  a  weak  form  of  the  even  parity  steady  state 
Boltzmann  equation  was  developed.  It  was  shown  that  because 
the  problem  equations  were  positive  definite  and  self-adjoint 
the  Rayleigh-Ri tz  variational  principle  and  the  Bubnov-Galerkin 
method  of  weighted  residuals  are  equivalent.  The  problem 
solution  was  formulated  by  using  a  trial  function  expansion  in 
bicubic  polynomial  splines  and  spherical  harmonics.  This  trial 
function  expansion  has  a  dual  basis  —  a  local  basis  in  space 
and  a  global  basis  in  angle. 

This  development  was  specialized  to  the  air-over-ground 
neutral  particle  transport  problem.  It  was  shown  that  a 
finite  element  space-angle  synthesis  solution  is  possible  and 
that  a  first  scatter  interpolation  source  can  be  used.  It 
was  also  shown  that  a  numerical  solution  can  be  achieved  and 
that  this  solution  technique  may  eliminate  rav  effects  and 
reduce  computational  costs. 

Recomnienda  L  ions 

The  preliminary  results  of  this  study  have  shown  that 
the  FESAS  method  can  produce  a  numerical  solution  to  the 
steauv  state  Boltzmann  equation  and  t  he  a  Lr-over-g round 
problem.  However,  because  of  the  time  constraints  on  this 


research  project  a  complete  development  and  evaluation  of  the 
FESAS  method  was  not  accomplished.  Therefore,  there  are  a 
number  of  recommendations  for  the  further  analysis  and 
evaluation  of  the  FESAS  method,  which  must  be  made.  Some 
of  these  recommendations  are: 

1.  To  enforce  the  boundary  conditions  at  the  air- 
over-ground  interface.  This  can  be  accomplished 
by  a  coalescing  or  stacking  of  the  nodes  (knots) 
of  the  bicubic  splines  and  by  using  a  Double- 

PM  approximation  at  this  interface. 

2.  Develop  or  refine  the  computer  algorithm  so  that 
a  comparative  study  can  be  made.  This  study 
should  include  a  comparison  of  the  computa t ional 
costs  and  accuracy  of  FESAS  to  those  of  Monte 
Carlo  and  discrete  ordinates.  Also,  a  determina¬ 
tion  should  be  made  as  to  whether  ray  effects 
have  been  eliminated. 

3.  Obtain,  if  possible,  a  closed  form  solution  to 
the  angle  integrals  in  Appendix  F. 

4.  Explore  other  ways  of  handling  the  discontinuity 
(at  .o  =  0)  of  the  space  integrals. 

5.  Use  other  spatial  trial  functions.  Lower  decree 
bi-quadratic  splines,  hermites  and  Lagrange 
polynomials  are  possible  candidates.  A  compari¬ 
son  of  the  results  winch  are  obtained  from  the 
use  of  various  trial  functions  can  then  be  made. 

6.  Examine  the  effects  of  an  improved  source  inter¬ 
polation  on  the  solution  accuracy  and  rate  of 
convergence.  An  improved  source  interpolation 
can  be  achieved  by  the  use  of  a  smaller  source 
(space)  grid  or  a  higher  degree  Lagrange  or 
Hermite  polynomial  interpolation  function 

7.  Extend  the  use  of  finite  element  space-ancle 
synthesis  to  the  solution  of  energy  dependent 
mult  : - r o u p  p r o b l e m s  . 
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Appendix  A 


Derivation  of  the  Scattering  Kernel,  Inverse  Collision 
Operators  and  the  Even  and  Odd  Parity  Collision  Cross  Sections 


Collision  Cross-Sections 

The  even  and  odd  parity  legendre  expansion  scattering  cross 
sections  o|  and  originate  in  the  derivation  of  the  second 
order  forms  of  the  Boltzman  equation  (Chapter  III).  A  defini¬ 
tion  of  these  quantities  has  been  given  elsewhere  (Ref  5:29). 
This  definition  is  repeated  below. 

Beginning  with  the  even  and  odd  parity  scattering  cross 
sections 


(A-l) 


(A-2) 


the  usual  computational  practice  (Ref  6:131)  is  to  expand  these 
scattering  cross  sections  in  Legendre  polynomials  as 

Z 

•3'LtL..  A  (A- A.')  (a -3) 

tel* 

or 

6s(rysi*SL) 


where 

o:’  ( r  ,  :i  •  H  ")  =  macroscopic  scattering  cross  section 

aS(r)  =  Legendre  macroscopic  cross-section 
expansion  coefficient  which  is  a 
function  of  position  (material) 

L  =  the  degree  Legendre  polynomial  expan¬ 
sion  which  is  used 


=/  6^(0- 2£+/  .  EfrZ-A') 


(A-4) 


P^ft-sV)  =  Legendre  polynomial  of  degree  Z. 

=  scattering  angle  (pQ  =  Cos6) 

Inserting  Eqs  (A-3)  and  (A-4)  into  (A-l)  and  (A-2) 


6. 


From  the  even  and  odd  properties  of  Legendre  polynomials 


17:223) 


even  function  if  Z  is  even 
odd  function  if  Z  is  odd 


therefore 


if  l  is  odd 


s  even 


and 


'2B6ZA')  if  1  is  cdd 


~£-o 


yjncA-A)  = 


o 


if  Z  is  even 


Eqs  (A-5)  and  (A-6)  can  now  be  rewritten  as 

Z 

=V 4iv  ■  B &■■*') 

I— j  4-  7T 


and 
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Z 


.jC'  -Z  / 


(A-5) 

(A-6) 

(Ref 

(A-7) 

(A-8) 

(A-9) 

(A-10) 


(A- 11) 


The  even  and  odd  parity  cross-sections  can  also  be 


expanded  in  Legendre  polynomials  as 


ina 


L _ i 

e*o 
z, 


(A-12) 


.  . . 


(A- 13) 

Comparing  Eqs  (A-”L3),  (A-12),  (A-ll)  and  (A-10)  it  is  apparent 
that 


and 


of(r)  =  \ 


OoC?)  = 


o  ?  f  r)  for  £.  even 
0  for  £  odd 


for  £  even 


<j^(r)  for  i  odd 


(A-14) 


(A-15) 


Therefore  of(r)  and  o®(r,.Q‘i1')  are  even  functions.  Also  c^(r) 


and  cr  (?,&•&')  are  odd  functions, 
s 


Scattering  Kernel 

In  the  development  of  the  even  and  odd  parity  forms  of  the 
anisotropic  steady  state  Boltzmann  equation  (Chapter  III),  it 
was  necessary  to  express  the  scattering  kernels  in  terms  of 
the  even  and  odd  parity  flux  (fluence)  components.  Following 
the  derivation  of  Wheaton  (Ref  5:11)  the  scattering  terms  can 
be  written  as 


where  the  integrations  are  carried  out  over  all  directions. 
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1 

1 


Because  a®  is  an  even  function,  meaning  f(x)  =  f(-x), 
it  follows  that 

em*.-*')  =  6/&A.A') 


C  A— 1 7 ) 


and  therefore 


With  the  even  parity  flux  defined  as 


Eq  (A-18)  becomes 


(A-18) 


(A-19) 


and  by  a  similar  derivation 


(A-20) 


(A-21) 


Inverse  Collision  Operators 

In  deriving  the  second  order  forms  of  the  Boltzmann 
equation  (Chapter  III)  the  and  Gu  operators  were  defined  as 

GWQ  -  &/&)  A )A>i')dyi '  (a-22 ) 

■vrr 

and 


( V 


Z.s  2 


1 


vrr 


(A- 2  3 ) 
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where  the  r  dependence  has  been  omitted  in  an  attempt  to 

simplify  the  notation. 

Usina  the  addition  theorem  (Ref  6:609) 


in  (A-12)  and  (A-13)  gives 


<£  3a-A')^  W 

A-o 


o  s*r-c  -<£^ 


Eqs  ( A — 2 2 )  and  (A-23)  can  now  be  written  as 


** L — 1  *  -%f?r 


(A-24) 


(A-25) 


(A-26) 


( A— 2  7 ) 


L 


syia) = &/a)  ^rTfj y/Ayj&^M 

The  inverse  operators  are  defined  as 


i  5  V-9W 

K -p 


K=[(T" 


where  it  is  mean 


t  that  if 


qYa)  =  R(A) 


then 


(A-28) 


(A-29) 


(A- 3d'., 


multiplying  Eq  (A-30)  by  /^Y^n(fi)dO  and  expanding  the 


operator  gives 


/&). y/yy&M, 


*^/W’  y/r 


I  />  'M 


(A-32) 


using  the  orthonormal  properties  of  spherical  harmonics  (Ref 


6:609) 


=  4C 

-hrr 


where  5  5 ,  is  the  Kronecker  delta  which  is  defined  by 


£.  = 


0  l  f  k 


1  £  =  k 


Therefore  Eq  (A-32)  can  now  be  written  as 


(A-33) 


(A-34) 


(A-33) 


*  4?  / 


(A-36) 


Rewriting  Eq  (36)  with  a  £m  spherical  harmonic  subscript 
and  substituting  into  the  expanded  form  of  Eq  (A~30)  gives 

R(A)  =  5^Xy/4?(A')yJaA,dA'  (A -r> 


Weak-Form  of  the  Functional 


The  functional  whose  Euler  equation  is  the  even  parity 
Boltzmann  equation  of  Chapter  III  is 


where  the  inner  product  <f,g>  is  defined  as 


and  *  means  the  complex  conjugate. 

The  minimizing  function  of  this  functional  is  the 
function  4*  which  is  a  solution  to  the  second  order  even 
parity  Boltzmann  equation  (Ref  10:169)  Therefore,  a  solu¬ 
tion  to  the  even  parity  Boltzmann  equation  can  be  found 
by  minimizing  Eq  (B-l). 

Another  more  useful  formulation  of  this  problem  is  the 
weak  form.  This  weak  form  can  be  found  by  imposing  the 
condition  that  a  function  which  satisfies  the  natural  boundary 
conditions  Eqs  (41)  and  (42)  and  the  even  parity  Boltzmann 
euqation  must  also  be  a  minimum  of  the  functional  (B-l). 

Let  the  functional,  (Eq  (B-L),  have  a  minimum  at  4'. 

Then,  for  all  n  and  c  where  o  can  be  arbitrarily  small  and 
ot  either  sign 


FM  <- 


X 


(B-J4 


where  '?  and  n  are  real  functions  that  satisfy  the  boundary 
conditions . 


Expanding  Eq  (B-l)  in  'i'  +  en  gives 


(B-4) 

(B-5) 

^j<A-  rz>  r !/))><& 

(B-6) 

r£J '<£■  #“£*■  VT/))  df 

(B-7) 

+ dr 

(B-8) 

+z^<¥,  d/r 

(  B— 9 ) 

(B-10) 
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(B- 11) 

-*£<*■  usuy  df 

(B-12) 

(B-13) 

'* £<<f,S3>  d? 

(B-l 4) 

sz{<^^sy^^ 

J«r 

( B— 1 5 ) 

^  <£  C //£■  ft  j  t/s 

(B-16) 

S2£Q  f  f‘/l.hj  y/  VcJstt/s 
-4-W 

(B-l 7) 

U- 

<S  '■»n' 

(B-  18 ) 

noting  that  (B-4)  +  (B~8)  +  (B~12)  +  (B-14) 

+  (B-io)  =  F  ( H' ) 

(B-19) 

6  b 


and  that  K®,KU,G^  and  Gu  are  self  adjoint  operators  (Ref 
10:174),  where  if  L  is  self  adjoint  then 

($,$*))  =  <f(A),  1 30Q) 

r  (  B-2 0 ) 

*  means  the  complex  congugate.  Since  n  and  f  are  both 

real  functions  ( B— 5 )  =  ( B— 6 )  and  (B-9)  =  (B-10) .  Therefore, 


collecting  terms  and  simplifying 


F(¥+£TC)  =  F(¥)+2£jj  (<sf-r%?/r&.v0} 


With  both  Ku  and  G®  being  positive  definite  operators  then  the 

2  2  • 
e  term  of  ( B— 2 1 )  must  also  be  positive.  Note  that  t  is  al¬ 
ways  positive.  Now  in  order  to  ensure  that  F('i  +  en)  >  F(T) 
for  e  5*  0,  the  e  term  in  equation  (B-21)  must  be  positive  or 
zero.  But  e  can  be  of  either  sign  therefore  the  coefficient 
of  e  must  be  zero,  that  is 


=  o  (B. 


Eq  (B-22)  is  the  weak  or  Calerkin  form  of  the  second  order  even 
parity  Boltzmann  equation.  A  detailed  derivation  of  equation 
(B-22)  can  be  found  elsewhere  (Rer  5:57). 


Appendix  C 


Derivation  of  the  Weak  Form  From  the  Galerkin 
_ Method  of  Weighted  Residuals _ 


It  can  be  shown  that  solutions  to  the  second  order  forms 
of  the  Boltzmann  equation,  by  using  a  variational  principle  or 
the  method  of  weighted  residuals  are  equivalent.  A  proof  of 
this  equivalency  for  the  even  parity  equation  is  outlined  below. 
However,  a  similar  proof  can  also  be  extended  to  the  odd 
parity  second  order  equation. 

The  starting  point  of  this  proof  is  the  even  parity 
Boltzmann  equation 

-A-rftriJi-vycr,*)  r  tpif) A)  (c-n 

and  vacuum  boundary  conditions 


(  C— 2 ) 
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r L '  K  (^)iS  (KrO  "t-  rp&s, 


■*o 


(C-3) 


for  xL.  ■  fv  O 


In  the  following  equations  the  r  and  9.  dependences  will  be 
omitted . 

If  a  trial  solution  T  is  assumed  where  ¥  is  a  linear 
combination  of  functions  such  that 

/V 


K 


(C-  4) 


then  the  Galerkin  method  of  weighted  residuals  requires 


)-<*■  vk  “S“X>  df  =  (<*  ‘5?*-  *  <1? 

-£((«  n) ds 

Substituting  Eq  (C~9)  and  (C-10)  into  ( C— 6 )  gives 

jjk-VTljfrk  wf}  ^<(63%7iyjdb 

f  <£ /^f/f  h)K“S“%. -  <(A  »)f<  “*■ 7>J  * 


(C-10) 


term  A  of  Eq  ( C— 1 1 )  can  be  rearranged  into 


(C-ll) 


IZ-A-V 


<5 


^ J  %d!n.<L£ 


( C  —  L  2 ) 


Using  the  boundary  condition  of  Eq  (C~2)  and  (C~3)  ,  Eq  ( C  —  1 2 ) 


can  be  written  as 


(/t.k)  (ftyd'k-  ( (kL.t 


d£ 


( C  —  1 3 ) 


AS  CAa,/>j 


'SL  ft<0 


Eq  (C-12)  can  also  be  written  as 


■<£  (//2-fy/ fydAds 
Jsl'T 


( C- 14) 


where  |j  moans  the  absolute  value.  Substituting  Eq  (C-14) 
into  (C-ll)  gives 

({/A.  ,*•<£?■ V,  r<£(/i.iii>.>7J/ia’s 

L  •  '  J  J.  /  _  '  ‘  - 


^sk-rr 


(C-13) 


Appendix  D 


Expansion  Properties  of  Spherical  Harmonics 


In  Chapters  III  and  IV  the  angular  dependence  of  the 

trial  functions  and  cross  sections  was  expanded  in  spherical 
harmonics.  Because  of  the  two  dimensional  angular 
dependence  in  y  and  xs  and  the  requirement  that  the  expan¬ 
sion  functions  should  form  a  complete  set,  the  expansion 
is  presented  with  m  and  i  subscripts  as  follows  (Ref  6:608) 


fa)  v  Y 


L _ ,  L — ' 

4?'  Z>  /xn  -r  —  £ 
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^  Y  Y  X 


(  D— 2 ) 


where 


~  J  J  xj  V x. )  d x 


(0-3) 


t?'7T7 


C  Hd 


(0-4) 


y>,  = 


X  'j  ^  \ 


(0-3) 


V**  YYuZ)/ 


(0-6) 


If  1 ( p , X )  is  even  in  the  angle  X  then  the  expansion  must  also 
be  even  in  and  therefore,  (D-4)  can  be  rewritten  as 


^  Yi,*  = 


(D-7) 


where  the  odd  iSinx  term  is  omitted.  Also  from  Eq  (D~5) 


(D-8) 


Therefore  for  an  even  expansion  Eq  (D-3)  can  be  written  as 


Jb 


(0-9) 


Jo  -/ 


^ ;/  /Up OIL zC-y  Tj  ^  ( [)- 10) 


For  uii  :ven  x  expansion  of  f(u,\)  Eq  (0-2)  can  therefore  be 


written  as 


L  f  r-  JL 

A/7KQt. 


(h-ro 
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L  JL 
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,'itiTi'  the  1  actor  of  .wo  has  been  inctuucu  in  Lhe  i  ,  ox- 

» :T> 

pans  ion  eoei f ioients  and  Lq  (e-fi)  forms  a  complete  sot. 

A  si.aii..'.  approach  can  bo  used  tor  the  cross-sect  io:i 


expans l on 


6UA-) 


<4  ^ 


A  /\ 

where  P^(P‘P.')  is  given  by  Eq  (A~24) . 
Eq  (D-13)  can  be  expanded  to  give 


£y'c^ 


and  from  Eq  (D-5) 


l  y.S')k  a)  r  kc  r)  lrS> 


Therefore  (D-14)  becomes 
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.;.  -i’zo  4r/*=  £> 

where  m  means  that  all  terms  with  a  m  =  0  subscript  must 


be  divided  by  two. 

Since  the  scattering  cross-sections  are  real  the  expan 
sion  of  Eq  (D-16)  must  also  be  real  and  with 

r  ] 

/'o~)  -  cT  V('*r)  •*"-  r-  Z. "  ' 

c  -J**  - 


—  /  ,  .*zl 


•iw  y- 


Eq  (D-16)  can  be  written  as 


1 


Note  that 


n(z-x') ^ni’Cxs7n%/-r  Sok,  D_19 


The  angles  y  and  \  are  shown  in  Fig.  2. 

Similar  derivations  would  give 
X  JL 


l-6i***6  ( H-21 ) 


By  inserting  Eq  (D-20)  into  Eq  (A-22)  the  even  parity  collision 
operator  can  now  bo  written  as 


an.!  from  Eqs  (P-13)  and.  (P-17)  t!io  inverse  odd  parity  collLson 
operator  Kq  (A-aO)  ,  can  be  written  as 


(D-25) 
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Appendix  E 

The  Synthesized  Boltzmann  Equation 

In  order  to  formulate  a  numerical  solution  to  the 
air-over-ground  problem  it  is  necessary  to  expand  the 
weak  form  of  the  even  parity  Boltzmann  equation  in  a  set  of 
trial  and  test  functions.  The  finite  element  space-angle 
synthesis  trial  and  weight  functions  of  Chapter  IV  will 
be  used  in  this  expansion.  Because  this  is  a  very 
tedious  and  lengthy  derivation,  the  more  obvious  albegraic 
steps  will  be  omitted.  The  starting  point  of  this  formula¬ 
tion  is  the  synthesized  even  parity  Boltzmann  equation,  Eq 
(58),  of  Chapter  III. 


and  they  are  given  by 
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Figure  9.  Surface  Normal  And  Particle 
Velocity  Direction  Vectors. 


and 


'5£7r  &>,»'){ 


P o  m( ^ )  are  t^ie  associatei^  Legendre  functions 


The  directional  derivative  is  given  as 

A 

A 


where 
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for  the  horizontal  top  and  boLto 
surfaces,  and 
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G 


A 


i  or  the  vertical  (side";  of  the 
nroM  ora  c  v  1  Lndor  ■  ' 


If  fi. 


n  and  n  are  considered  to  be  unit  vectors  then 
r  z 


from  Fig.  7 


RB  -  EB  =  CF  =/ 

and  therefore  for  the  top  surface  of  the  cylinder 

/£. **  Rg*££>*£asQ  */*■ 

and  for  the  bottom  surface 


(E-10) 


(E-U) 


r  -  l/s 

also  CD  =  FB  =  projection  of  H  unto  the  x-y  plane 
and  therefore, 

CD  =  HB  =.  EB  SinO 


(E-12) 


=  -\j  1-Cos20  =  -yj  1-y2  (F-13) 

and 

/l  np=CV  *£? 

(E-14) 


In  Appendix  D  the  even  and  odd  operators  were  given  as 


and 
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7D 


where  m  means  that  all  terms  with  a  m  =  0  index  must  be 
divided  by  two,  and 


Also  the  inner  product  is  defined  as 


(  E— 17) 


d ^  forrealf  ( E— 1 S ) 

Using  Eqs  (E-16),  (E-15),  (E-14),  (E-ll),  (E-12)  and  (E-6)  and 
noting  that  T  r  .  Q  „  ,  B-  and  B-  are  all  real  functions,  F.q  ( E  —  1 ) 
can  be  expanded  to  give 
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(  E— 35 ) 


(E-36) 


( E— 3  7 ) 

(E -38) 


=  C i'' Sl/A 
*  /*  2)/=>  J  2*» 


v'E-40) 


AfgA-Sa'^' 

'ClTf 


(E-41) 


where 


dC 
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(E-47) 


and 


Eqs  (E-19)  to  (E~36)  is  an  expansion  of  the  first  term  of 
Eq  (E-l).  Eqs  ( E—  3 7 J  and  (E~38)  is  an  expansion  of  the  second 
term.  Eq  (E-39)  is  an  expansion  of  the  third  term.  Eqs 
(E-40)  to  (E-46)  is  an  expansion  of  the  right  hand  side  of 
Eq  (E-l).  o^,  ot  and  o^r  are  functions  of  z  and  they  must 

be  included  in  the  spatial  or  dv  integrals.  In  cylindrical 
geometry  with  azimuthal  symmetry 


A\/  -  2  Tf/'d/'c/  z 


(E-49) 


and  fds  means  an  integration  over  the  surface  of  the  problem 
cylinder,  Fig.  9. 

For  the  air-over-ground  problem  with  an  exponentially 


*3 


varying  air  density 


6&>  =4^><r  *** 

(E-50) 

(E-51) 
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where  ot(o) ,  o^(o) ,  o^!  o)  are  cross-sections 

of  air  at  sea 

level . 

Z  =  the  height  above  sea-level 
sh  =  atmospheric  scale  height  ~  7km 
In  Eqs  (E-19)  to  (E-46)  the  integrals  are  separated 
in  the  space  and  angle  variables.  These  are  double  integrals 
in  space  and  angle.  However,  they  can  be  separated  into 
single  integrals  of  the  p,  x>  P  and  z  variables. 


Appendix  F 


Angle  Integrals  of  the  Synthesized  Second 
_ Order  Boltzmann  Equation _ 


An  expansion  of  the  even  parity  Boltzmann  equation  has 

produced  twenty-eight  integral  terms  (Appendix  E) .  By  a 
further  expansion  and  separation  of  the  integration  variables 
twenty  distinct  single  angle  integrals  are  formed.  These 
angle  integrals  are  dependent  on  the  degree  of  the  spherical 
harmonic  trial  function  expansion  and  independent  of  the 
problem  parameters.  They  can  be  evaluated  once  and  thereafter 
used  as  a  part  of  the  problem  input  data.  In  this  research 
project  these  integrals  were  numerically  integrated  for  each 
combination  of  the  m,  k  and  n  expansion  subscripts.  They 
were  then  stored  as  a  matrix,  and  selected  products  were  used 
to  produce  each  of  the  twenty-eight  angle  integrals  of  Eqs 
(E-19)  to  (E-46)  in  Appendix  E. 

These  twenty  integrals  are 
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A  numerical  evaluation  of  all  twenty  integrals  was 
carried  out  for  a  third  degree  spherical  harmonic  expansion. 
This  evaluation  showed  that  these  integrals  are  equal  to  zero 
for  many  combinations  of  the  Jim  and  kn  subscripts. 

Integrals  (F-10)  to  (F-14)  are  a  part  of  the  surface  inte~ 
gral  term  which  has  been  partitioned  into  the  outward 

A  «*S 

(fi*n  >  0)  and  inward  (ft*n  <  0)  directions.  This  partitioning 

was  incorporated  into  the  weak  form  derivation  of  Appendix  C. 
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Appendix  G 


Space  Integrals  (Bicubic  Splines)  of  the  Synthesized 
Second  Order  Even  Parity  Boltzmann  Equation 


A  trial  function  expansion  of  the  spatial  flux  dependence, 
in  the  weak  form  of  the  even  parity  Boltzmann  equation,  has 
been  carried  out  in  Appendix  E.  Bicubic  polynomial  splines  in 
the  P  and  z  variables  were  used  to  form  a  tensor  product 
space.  These  splines  are  twice  continuously  differentiable  and 
have  non-zero  integrals  /^B^(x)B j (x)dx  for  all  |i-j|*4. 

After  a  separation  of  the  p  and  z  variables  of  integration 
seventeen  distinct  integral  forms  are  produced.  These  integrals, 
which  include  the  source  integrals,  must  be  evaluated  over  the 
entire  problem  domain.  The  space  integrals  of  Appendix  E  are 
selected  products  of  the  following  seventeen  single  integrals. 


v 

(* SM  Jr 


(G-l) 


(G-2) 


/  3.  G*)3.  C/°)  A/? 


a?  JA  /A 


(G-3) 


CP 

Jr 

l 


(A  .  Kcr) 

Jr  //- 


(G-4) 


(G-5) 


'  j>-  "*■ 


(G-6) 


(G-7) 
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Z/Z)Jz- 


&,{*)■  e~z/sA  dz 


o 

rH 


JL 


Zfz).&6Cz)  dz. 


The  Source  Integrals 
'X 


(  HjOtes 

V 

(  Hjfrf  A  p 


& 


&(P)  ft{A>)  /?J/° 
&r  J 


B  <z)  H  Cz)  ^ 
J*  / 


Jz 


H 


U*Jz)J  »/*>  *  ZA"Az 

O-Z- 


hi 


where 


B(z) 

B(x) 

sh 

R 

H 


3/z)  ■  H.  (z)  Jz 

J*  ' 


cubic  polynomial  z-spline 
cubic  polynomial  p-spline 
atmospheric  scale  height 
outer  radius  of  the  problem  cylinder 
problem  cylinder  height 


(G-8) 

(G-9) 

(G-10) 

(G-ll) 

(G-12) 

(G-13) 

(G-14) 

(G-15) 

(G-16) 

C  G— 17) 
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H(z)  and  H(p)  are  the  source  interpolating  functions  (linear 
Lagrange  polynomials). 
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An  Expansion  of  the  First  Scatter 
Source  in  Legendre  Polynomials 


In  Chapter  IV  the  first  scatter  source  was  defined  as 

■S(P,z,A)  --  (H-i) 

where 

($i ^ C r , z  ,Q  ')  =  direct  fluence  of  Chapter  IV,  Eq  (74) 
os(z,r>.fi')  =  scattering  cross-section 

The  usual  Legendre  polynomial  cross-section  expansion  will 
now  be  carried  out.  Also  the  even  and  odd  parity  first  scatter 
source  expressions  of  Chapter  IV  will  be  derived.  Expanding 
as  in  Legtndre  polynomials  and  using  the  addition  theorem 
(see  Appendix  D) 


z  -e 


C -O 


where  m  means  that  all  terms  with  a  m  =  0  subscript  must 
be  divided  by  two,  and 


6*a)  =  Sjcd)  e 


-  */sA 


From  Fig.  5  and  Fig.  2  it  is  apparent  that  y  ■*  =  yd 


and  x  *  =  0 
Therefore 


S(A  2t  A)  ,  <2  <pj  tfz,A')  £  -  ** 

*  5~w;£^s'~  * 


(H-4) 


£■=<>  ****<? 

and  using  the  identity  (Ref  18:96) 


(H-5) 


and  a  little  algebra,  the  even  and  odd  parity  first  scatter 
sources  can  be  written  as 


j. 

2. 

=  & 


'SC'S*,*)  -SCte-ti 


"42  ~C>  ***£> 


j  (H-6) 


and 


f  S(^zy  ~/£)\ 


*-) 


e-  M 


■^-e>  4*  *. 


= 


7^)  ~P<C<*cZ)  X) 


-erM 


(  H  —  7  ) 


cy> 


A  Derivation  of  the  Total  Particle  Fluence 


In  Chapter  IV  a  trial  function  expansion  of  the  even 
parity  angular  particle  fluence  was  given  as 

fff  fkfrPQ&Q 

L - 1  L - 1  L _ i  C_ _ i  ' 


/>=/  •V'r/ 


0  /frt-0 


where 


(i-D 


(1-2) 


4  ^ 


4^- -v  \ 

V  ^  ^  <^V  7*)/  (1-3) 

The  A-  .  „  mixing  coefficients  are  obtained  from  a  numerical 
i  ,  j , £ ,m  ° 

solution  of  the  second  order  synthesized  Boltzmann  equation  of 


Chapter  IV. 


The  angular  even  parity  fluence  is  also  defined  as 


^4^ yU^')  v-  >v  (i-^ 


An  integration  of  Eq  (I~4)  over  all  directions  gives  the 
total  even  parity  particle  fluence  T(p,z),  and  also  the  total 

particle  fluence  ;(e,z't.  This  is  because 


[W )dt t  -  z ) 

/dir 


d-3) 


Therefore 


Eq  (Z-I)  will  now  be  integrated  to  give  the  total  particle 

£1uence  at  position  (P,a).  Using  the  orthogonal  properties 
of  Legendre  polynomials  this  integration  is  carried 

follows. 

The  aero  order  associated  Legendre  function  rs  deftned 


as 


'Oj0 


(1-7) 


.  .  .  Vn  rr_n  5v  Fq  (1-7)  and  integrating  over  all  o  and 

Mul tiolvmc  fq  (I  D  D>  KL  ' 


94 


where 


'2ir 


=1 


and  (see  Appendix  D) 


&  for  m  f  0 

•2  ^  for  m  =  0 


(1-10) 


-  2£  ~  -  2 //err* 


(i-ll) 


therefore  Eq  (1-9)  becomes 


2  7T.2//^~rr  -J^TT 


(1-12) 


and  the  total  particle  fluence  is 


JT2 


where  JL^  -?n  -  O 

The  angular  particle  fluence  is  given  by 


2*, *)  +  X(/>Z,  X)  (i-u) 


where  x(P»z»»0  i-s  the  odd  parity  fluence,  which,  is  defined 
in  terms  of  the  even  parity  fluence  and  source  by  the  follow¬ 
ing  expression 


V/>7/f)^  /r  V> 


(1-15) 


1 
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Appendix  J 

Computer  Subroutines 


A  computer  program  has  been  written  to  solve  the  synthe¬ 
sized  Boltzmann  equation  of  Chapter  IV.  This  program  is 
designed  to  use  a  trial  function  expansion  in  bicubic  splines 
and  spherical  harmonics,  and  to  perform  a  first  scatter  source 
interpolation  using  linear  Lagrange  polynomials.  The  program 
in  an  assemblage  of  several  subroutines  which  collectively 
perform  the  following  tasks. 

1.  Computes  all  single  space  and  anele  integrals. 

2.  Combines  the  single  u  and  \  angle  integrals  for 
all  combinations  of  the  spherical  harmonic 
expansion  subscripts. 

3.  Combines  the  single  o  and  z  integrals  for  all 
combinations  of  the  bicubic  spline  expansion 
subscripts . 

4.  Assembles  the  coefficient  problem  matrix. 

5.  Computes  and  interpolates  the  first  scatter 
source. 

6.  Assembles  the  source  vector. 

7.  Checks  for  symmetry  and  diagonal  dominance. 

8.  Solves  for  the  A ; , ]  ;  n  expansion  (mixing) 
coefficients  by  the 'method  of  successive 
over-relaxation. 

9.  Solves  for  the  total  particle  fluence. 

A  ten  point  Newton-Cotes  single  integration  routine  was 
used  to  numerical  lv  integral  •'  tin  thirty-seven  integrals 
of  Appendices  F  and  G.  l'his  integration  routine  is  an  in-hou 
subroutine  of  the  Air  Force  Aeronautical  Systems  Division , 

V.'ri  .-hi  -L’utterson  Air  Force  Base .  The  overall  program  logic 
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has  been  written  in  a  manner  whereby  this  integration  routine 
could  be  used.  The  program  is  written  in  Fortran  V. 


Listing  of  Problem  Subroutines 


PROGRAM  MR  IN 

£|  I  MENS  I QN  CL  <  0 :  3  -  0 :  3 '•  *  S H  •  £  >  ,  DUS  H  ■  0 :  3 *  0 : 

1  -z'  *  Z*  *  L  H  V  ->  9  'z*  J  -Z*  I.1  *'  •  •  '*£*  *  *7  1  •  .  L  V  j  ' 

D I  MENS  I  ON  PNZ  •  -  i :  7:.-  » PNP  •  -£ :  7>  *  FNX  ■-£:  7>  *  S  PZ  '7.7. 
17»7»6>  » KSEVL  C  0 :  4>  *  SVZ  7  .  5  *  3>  *  SVX  '7.5.  3>  «  S  VP  ■:  7  *  5 « 
1  CO  .  AS  <50>  SCO  *  XSOBL  *•' '0:  4>  •  XS'  CO;  4> 

COMMON.  R 1  'F'NK  .  I  A£  RS  H  ■  I  TYF'E  <  IT'S*  I  TTY 
COMMON  ."BE-  "ISH  E3-  IP 


. PL '3* 3* S> . DU PL  1 

*  3 1 P  P  1  7  >  7  *  €• . 1  *  s  F1  ’ 1 

*  RPD  '.'5*  5*  -SO  *  SB  •  5 


♦  THIS  ROUTINE  SOLVES  THE  EVEN  PARITY  ANISOTROPIC  BOLTZMANN  ♦ 

♦  EQUATION  AND  THE  AIF'-OVEP  GROUND  PROBLEM  FOP  THE  REACTION  ♦ 

♦  PATE  FLUENCE  BY  CALLING  A  NUMBER  DF  SUBROUTINES  AND  USING  A  ♦ 

♦  TEN  POINT  MENTON  COTES  INTEGRATION  ROUTINE. 


READ  1  ♦  *  * .  END  =  1  £  01  DZ  *  Z £  .  Pel .  DR  *  Z I  *  P  I  *  AsH  *  L MAX 
READ  c  ♦  *  ♦  *  END=  1  c'  u  1  YD  *  HB*  S I GT « t«l 
PRINTS*  PROBLEM  INPUT  DATA 
PRINT** 

PRINT*.  DZ=  .  DZ*  DF*  *  DR."  P£=".P£*  Z£=  ■".!£-  ASH* 

1 .  AS  H .  Z 1  =  *  Z 1  .  P 1  =  .  R 1 «  F'ELAXAT  I  ON  F  Ac  TDF:  = '  .  M 

PRINT*. 

PRINT*.  LEGENDRE  EXPANSION  P  •  LMAX>  = *  LMAX 
READ • ♦.  *> END* ISO  ■  CKSEVL • I > * I=0*LMAX> 

READ  END*  1  £ 0 >  CX S'ODL  •  I  1  *  1  =  0.  LMAX> 

READ  >.♦.♦*  END*  1  £ 0>  CX I  1  I > . 1  =  0. LMAX' 

PRINT*. "YD*" . YD*  "  HB*  .HB.  SIGT='.SIGT 

PRINT*."'  ■ 

PRINT*.  LEGENDRE  ODD  EXPANSION  SCATTERING  CROSS  SECTION.  ALL  ODD 
1  MOMENTS  ARE  ZERO. 

Pr I  NT* . 

-C  i  '  «  '  ’  ~  E!  ■ ! .  I  •  *  J  ~  :  i  »  {__  t’'  H '  * 

*-■  I  p{  T  ♦  * 

PRINT*.  LEGENDF  E  ODD  PARITY  EXPAN.IOM  SCATTER  IN'?  CPOSS-IECTICN. 

ILL  -L'  EVEN  MOMENTS  APE  ZERO. 

PRINT*. 

PRINT*.  '  .  XSODL  •  I  •  .  1  =  0.  LMA:  A 

PRINT*. 

F  F I M  T  * .  LEGENDRE  EXPANSION  CROSS  SECTION  'ODD  AND  EVEN', 
pc  rrir., 

r  =  T  r  ?  T  ♦  , 

j-  p  t  r  j  •  • 

7  -  •  ‘-’K  t-  l  ■  ♦  •  LP^!  -  s  -l' 


y 


y  n  n 
rszT  , 


9S 


J 


PRINT*.  THE  SPHERICAL  HARMONIC  C QEFF I C I  ENTS  1 L  «  N 1  . 

PRINT*. 

no  5  l=  u  -  lmax 

PRINT*.  ■  .  C L  >: L  *  N >  .  N  =  0 .  L 

PRINT*. 

PRINT*. 

K  1  =  0. 

X£=8*ATAN  >1 .  :> 


TOL=. 0000 l 


FOUNT =£00 
ITYPE=1 

PRINT*.  ANGULAR  INTEGRALS  OF  SPHERICAL  HARMONICS  FOP  THE 
PRINT*.  GALEFFIN  SOLUTION  TO  THE  EVEN-PAR  I  TV  FORM  OF  THE 
PRINT*.  •'  NEUTRON  TRANSPORT  EQUATION. 


PRINT*. ■  THESE  APE  INTEGRALS  OVER  2*PHI  DF  THE  PEAL  AND  OF 
PRINT*.  "  PART  OF  EXP  • -IMX -  .  EXPANDED  IN  SINES  AND  COSINES. 
FRINT*. 

DO  15  I SH= 1.1c 


L  ALL  SF'HER  1  LMAX .  _  H  1  U .  LN  I  S  H  1  .  X 1  >  Xc- .  TDL  >  F  OIJNT.J 
PRINT*.  "  ANNULAR'  ‘HARMON I L  1  INTEGRAL3"  .  ISH 


DO  10  M= 0 . LMAX 

PRINT*.  1  .  _  H 1 M . U «  I S H 1 . N  =  U . L  M A  ‘  . 1 

PRINT*. 

PRINT*. 


CONTINUE 
Xl  =  -1 .  0 
X£=  1 .  0 
ITVPE=£ 

PRINT*.  THE  ASSOCIATED  LENGENDRE  POLYNOMIALS  INTEGRALS  FOP 
PRINT*.  GALEFFIN  SOLUTION. 

PRINT*.  THESE  APE  THE  POLYNOMIALS  F'L  -  L  .  M :-  -ASSOC.  LEGENDRE- 
PRINT*.  ■  INTEGRATED  OVER  THE  INTERVAL  -.-l.+l>. 

PRINT*. 

DO  £5  I P  =  1 .  '£ 

CALL  POLY  - LMAX. PL  -  1 .  1 .  IP> . X 1 . X£. TOL . FOUNT  - 

PRINT*.  ASSOCIATED  LEGENDRE  POLYNOMIAL  INTEGRAL3 '' .  IP 

DO  £0  1 1.1=1.  FT 

=-p  TNT  *  .  <  FL  T'-N  IT  -  IF  -  •  IT  =  1  .  F  T  - 

PRINT** 

f '  P  T  -  T  ♦  . 


CONTINUE 

CALL  iNuLE  -  LMA!  N  L  L  *  ,h«  PL «  L  A .  S  I  NT  »  XS  EVL .  Xs  DDL  .FT.  X  V  '- 
ITT',  =1 

NX= - Z£-Z 1  -  DO* 1.1 
II  07  1 1  I  =  f 1 .  f  F  + 

C  fr  '  .  t  ,  =  T  !  +  T  *  r-S 

-ns  !.■  =Rff' .  r  • 

■  -  '  S  -  -  - '  r  - 
[ TiP F  = ! 


EVEN 


THE 


e  f  r-  r  ».  . .  m  j  i  1  cr  ~  ~rp  'urc  ;  .  .tegc-p-L 


—  r  i  r  r  >.• 
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iti  =  i 

print*.  mhtpi::  vrluei  fop  z-iprce  zoupce  integprli 

PRINT*. 

DO  OP  1=1 <3 

c hll  : dupce 1 i vz . nzs - tol - f  ount ■ 

CONTINUE 

n:  :=  •PE-Pi  '  df+i .  i 

DO  6  0  I  =  0-  fi'  .'+c 
PNM ' I ■ =P 1 + I ♦DP 

if- pn: : i -< . ec.  o.  o > fn: ■  •  i  ■  =  i  .  oe-p 

PNP  I  >  =  F'NK  I 

np:  =n::+p 

NT  =  1  T  +  F  T ♦  '  •  NZ  Z  -  1  '  +NZI  ♦  ■  NFS-  1 
PP I NT*. 

PFINt*. 

PP  I  NT ♦ *  THE  TOTRL  NUMBER'  OF  TPIRL  FUNCTION'  =  .NT 
PP I  NT* .  THE  NUNPEP  OF  LEGEND PE  FUNCTION!  =  HT 
PRINT** 

PRINT** 

PRINT**  PROBLEM  GEOMETRY  INPUT  DRTR. 

PRINT**  NODE  1  Z*  P  '  MEIH  C  DOP'D  I NRTE  '  Z  *  P  • 

PRINT** 

DO  :-5  NP=0*NZ'I-3 
DO  6 O  ML  =  n*NP.-;: 

F'  P  I N  T  ♦  »  1  *  M P1  *  *  *  Ml  *  1  *  *  F1  N Z  1  M F  1  *  *  p N P  1  M C 

print*. 

PRINT** 

I  TYFE=c 
it:=o 

PRINT**  MR TP  I ' :  VRLLEI  FOP  P-IPRCE  INTEGFRL I . 

PRINT*. 

DO  1 1  1  =  1  *t 

CRLL  :  'FLINE  •  IFF.NF  ;  *  TOL*  FOUNT  • 

CONTINUE 

it:=z 

PRINT*.  MRTF-i::  VRLUEI  FOP  P-ZF'RCE  Z OUPCE  I NTEGF'RL'I  . 
PRINT*. 

T<  p  <  I*.  T  —  1  “ 

C  R1  -  OF  P  :  •  .YR  *  HR  I  *  TOL  -  •  flUNT'* 

cr.r«Tir  sue 

TO  COMPUTE  THE  INTER PCLRT ION  POINT  VRLUEI 
L  RL L  D  I  F  E L  T  1  L.MR!  *  RR  D  *  F'tL.  *  F  NP'  *  NF’  -  *  NZ  .  *  .  I G  T  *  H F  *  YD  *  H  _  H  7 
TO  COMPUTE  THE  PROBLEM  I OUPCE  VECTOP 
RL  L.  E'VE L  1  RP D »  _  *  .  VP'  *  _  £  *  NZ  _  *  NFy  *  LMR.  *  L  R  *  _  1 1. T  *  1  T  *  N  T  1 

TO  COMPUTE  THE  ITIFFHEII  RTP I . 

1  RL  1-  .  T  I FR  « i_  h  *  _  F  P'  *  _  F' _  *  T L . *  NP'  _  *  I GT  *  H  I  *  F'NF  *  F'MZ  *  L MR Y  .  N  T  1 

Y  ;hE'!  'UR  IT  IFMEI  I  MR7  R  T ' 

-,l_  I_  ■  up,-  V  ,  h  -  ,  f-4  T  ' 

TO  I  OL’-E  THE  PROBLEM  MRrP  10. 

‘  r.L  'R  .  '  r  ,  t-  ‘  .  ■  ■  T  . '  " 

rr;  '  ON*1 i  •  rE  thp  _  ~  *—  ;  -L  rlRR  Fl.uy  gt  th£  f<CDE  '  . 

F  P  I  N  T  ♦  .  i  HE  rzr.iL  r  Ri.  I  RL  hi  /  EPhi.  h  p  D  I  P  E  l  1  R  L  CL  . 

rr  c  [  *  i  T  *  ,  ‘r.Cii'E  ,;C!NT  ~  ■  r  .  tht-g  •rpTf'Cft. 

L  ,*•  I  c  c  ■  T  rr  r-  ■  •  r  r  c  IT,-  T  r  -  r  hTTc  r  , 
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m^=  •  nz_-c ■  ♦c 

MP= ' NPC-C' ♦£ 

DZ  = ' ZE-Z1 ■  HZ 
Dr  -  ■  PE-P1'  HP 
DO  110  1=0. MZ 
DO  110  J= 0 . HP 
P=P1+J»DP 

IF'P.LT. PNP C 0 > > P=PNP >  0 • 

Z=Z 1 + I ♦DZ 

l  Fill  tfluy  <.  pop  ,  fnz .  xe  .  p  .  z .  i  t  ,  nz  :.  ,  np  _  ,  tf'H i  ,  dfh  i  >  he «  v i gt .  -  h .  yd  1 

TFH I =DPH I +TPHI 
PP I  NT ♦ ,  ' 

PP  I  NT ♦  Z ...  F' '  ,  TPH I  *  >  DPH  I  * 

1  -TFH I 

CONTINUE 
GO  TO  130 

PP I NT ♦ ,  PPEMPTUPE  PPOELEM  DRTR  END 

C  TOP 

END 

:  U  DP  OUT  I  HE  I PHEF  •  LHRY ,  DUE  H ,  X 1 ,  X£ .  TDL .  k  CUNT  :• 

dihen: ion  du:h ■ o: 3. 0: c ■ 

C  OH  MON  E 1  M  -  N  ■  EE  I  C  H  I  TYPE ,  ITS.  I  TTY 

♦  thi:  pdutine  nunepichlly  integrate!  the  cine:  end  cocinec  .tingle 

♦  INTEGRAL.!  OF  THE  RNIIOTPOPIC  EVEN  PARITY  BOLTZMANN  EOLATION  I"  * 

♦  CUE  OF  R  TEN  POINT  NEMTQN-C GTE!  INTEGPRT ION  PCUTINE.  ♦ 


IF ■  I IH. LE. ?• GO  TO  E 0 
IF  1  I  IH. EG.  10>  THEN 

:  1  =  0. 

:  :e=e.  ♦rtrn  ■  1 .  • 

GO  TO  EG 

el:e  if  'Ich.eu.11'  then 
:  :z=e.  ♦rtrn  .  1 .  . 

: :  1  =f#rt  mn  'j.'.' 

GO  TO  EG 

ELCE  IF  'ICH.EO. IE'  THEN 
: :  1  =r.#RT  rn  '1.' 

,E  =  -  *“  ■  !  . 

END  I  F 

TINT  I  "I  £ 

DO  ?M  H=0.  l hr:  : 

DO  CO  N  =  0  •  LHRY 

CRLL  '''HAD  1  ."1  . :  :£.  TOL .  DUCH  •  M.  N  •  .  I  OUMT) 
I  F  •  r»J  :  H  '  H .  N  ■  .  L  T  .  1  .  me-  l  0  '  DU  I H  •  M .  N  •  =  0  . 

continue 


*  i  .t-cgt  r  fNc  r  O!  Y  ■  L ' n''  •  I !  1  . ' .  I . :  :E  •  TIT  •  >  C'  ;N 

d  I mfh  ; :  rzr<  ri ,r  .  ,  ■ 

zMvrT  :  :  r'!.!'  i  _•  ir.  =e-  it:,  it 


»  riJ  I  r.  •"!!_:  r  t  ‘  ,r  .  ; ’.-it  i  |  r  e  "  tji  Z1-  TJM  -*  ~  r"  JRTEE  LF‘;£' 

♦  POLYNOMIAL  I  I.IHICH  Hpp  a  cret  of  THE  c.'NTHF:  I  zed  fvfn  prptt- 


■  ;r  r  {  z*  p  •  *  j 


::  ICM 
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IF ■ IF . LE. E ■ GO  TD  ?n 
IF  'IF'.  EO .  ?:>  THEM 
:  :i  =  o. 

:■  :g=  l .  o 
GO  TO  30 

ELCE  IF  'IP.  EI'.  3 '  them 
■ l  =  o .  o 
:  :e=-i  .  o 

END  IF 
CONTINUE 
1 1,1=0 

do  4  0  f  =  o .  lmh.l 
DC  40  M=0»I 
I T  =  0 
11.1=11.1+1 

DO  40  l  =  0, lmh:: 

DO  40  M=0,L 
I T= I T+ 1 

CALL  C'MhD  'Ll.'  3  .  TCL  .  DUEL  •  II.U  IT'  .  1  DUMT  • 

IF  ■  D'JPL.  •  Ii.M  IT'  .  LT.  1 .  0E-  1  0>  DUF'L  •  I M  *  IT  '<  =  0 . 

COMTIMlJE 
c  ETUF'M 
EMD 

FIJMC  T  I  CM  V'"' 

D I  MEN  T  I  CM  PM.' '  ■  -  3:  T  ■ 

COMMON  El  M.M-E'E  IIH  Dj  IF  E4  L  >  I  I  TVF'E «  I  T  I  «  I  TTY 
COMMON  FI  PM'  .  I  He  Pl'H  h?  JP.  JT  H4  1T.IH.LF 


♦  THi:  FLMfTICM  i  :  I.LED  BY  THE  TEM  POINT  MEMTOM  COTE  '  I MTEGFFT I  Of 

•  c'CMT ime  -i'.upd .  to  ielect  the  ihtegfhticn  function:.  it  u:e:  the 

+  I  M c  0 F' MA  T  I C  M  IM  THE  Af.C'-E  COMMON  FLOC1  :  TC  DETERMINE  THE  FIJMC  TIC 
+  I.IHICH  I!  BEING  I MTEGF  FtTED . 


IP  ■  ITT’.'.  EC  .  1  'GO  TO  £c 
IF  ■  IT'.'FE.  EC.  E  '  GO  ^C  3  0 
IP  '  I T  n . EC.  1  ■  THEM 
7  =  f  d  *  1  n  *  .  f'iv  ■ .  *r  o  * 


,  - .  c  .  *  ■  c C  *  .  0  +  ■  -l'+  I L  • :  '  +  I  I M  1  m+:  •  ■  +L □ .  -  N*L 

FETIJFM 

el:e  ip  ■ i :h. pp. : . . i :h. ec . i 
v=c o :  ■ ' '  • +i.  c i ' m+'  .  ♦re :  ■  m+:  . 


0.  OF  .  I  :h.EO.  1  1 .  OF.  I  CH.  EC.  13  ■  THE!-1 


V=C QS  <X>  ♦COS'  *'Z  IN  1  M^X> 

RETURN 

EL  I E  IF  •  I SH.  EC'.  8  •  THEN 
Y=COS  ■.'N«>  ♦SIN  1  •' 

RETURN 

ELSE  IF  aSH.EQ.9-'  THEN 

Y=  'COS  OO  ♦COS  '  N^X  '■  -N^  S  I N  X>  ♦S  I N  ■  N^X  •  >  ♦!  I N  •:  M* 

END  IF 

RETURN 

CONTINUE 

H=i-:  :♦♦£ 

I F  '  h  .  L  E .  1 .  u  E  —  1  U  ft=  u . 

IF  ■ IP. EQ. 1  •  THEN 

Y-=H»F LF  i-;x»  K  •  N  ■  ♦PLF  CX*  L  <  MS' 

RETURN 

ELSE  IF  aP.EQ.£>  THEN 

Y— SOF' T  1  fl  ♦'•‘♦PLF  1  X ’  K  *  N  1  ♦PLF  1  X »  L  *  N  ’ 

RETURN 

ELSE  IF  < IF. EQ. 3  S'  THEN 
—  ■_  L'R T  1  H  1  ♦PLF  1  X »  P  ?  N  1  ♦PLF  1 X »  L  »  M  > 

RETURN 

ELSE  IF  'IP.  EQ.  4"*  THEN 
Y  =  !  '♦♦d*F‘LF  1  X  •  F  i  h  1  ♦PLF  ' . »  L »  H  -' 

RETURN 

ELSE  IF  ■  IP.  EQ.  5.  OF.  IF.  EQ.  7.  OR  .  IP.  EQ.  S>  THEN 
♦F  L F  1  P  ’  N  >  ♦PLF  1  X  •  L  »  M  1 
RETURN 

ELSE  IF  'IP.EQ.t'  THEN 

y-plf  ' f  .  to  ♦plf  - x«  l «  n  < 

END  IF 

RETURN 

CONTINUE 

IF  ■'  I  TYPE.  EQ.  S  '  GO  TO  S  O 
IF- ITS. EQ. 1 • 80  TL  ! EC 
IF  ■ I.EC.l'  THEN 
M=t 
N=  l 

i  -  n  t  n  p 

else  a  ixsr.s  ’hen 


n=s 

QC  TO  70 

ELSE  IF  a.EQ.S'  THEN 

n=r 

QD  TO  70 


'"0  th  —  r, 
'-■'DIF 
-  rnr  i  ni  t 


in  3 


i .  a 


IF  •  ITS.  EC'.  £■'  GO  TO  130 
IF  ■: I .  EG1.  1  >  THEM 

M  =  3 
M=3 

GO  TO  70 

el:e  if  •i.ec'.e*  them 

M=3 

i 

GO  TO  7  m 

el:e  if  •i.eq.3.-'  them 

M=  1 
M=  1 

G 0  TO  7  U 

EL  IE  IF  •  I  .EC'.  4  •  THEM 

M=3 

M=  1 

GO  TO  70 

EL :  E  IF  a.EQ.5>  THEM 

M=  1 

M=1 

b 0  TO  7  1 1 

EL  IE  IF  ■  I .  EQ.  6>  THEM 
M=  1 
M=  1 
END  IF 

iPi  =  :f-f < m.i  - f-m: : •  jf-g  ■  ? pm: : .  jp-£> . pm: •  ■  jp-i ■ . pm:-: •  jp> 

I F' _  F'  F  1  f  I «  1  —  J T  *  F' M : : '  J P'  +  J  T  —  1  »  F'  r  i : :  1  J F'  +  J  T  ~  c . 1  '  F' M : :  1  J P'  +  J T  —  1  '■  <  PM 7  '■  J  P  +  -l T 
1 ' 

IF  1  I  TYPE .  EG.  c'  ■  GO  TO  100 
IF  1  I . L E . 3  >  THEM 
V=E.VP  '  .'•••;  H  i  H  1  ♦  I P 1  ♦  ;•  F  3 
FETIJPM 

ELIE  I P  1  I . E  0 . 4 1  THEM 
v=e:  P'-':  ft: H  •  ♦  i P l  ♦  i P£ 

FETIJPM 

ELIE  IF  1  I . EO .  j 1  THEM 
V=  I P  1  ♦  '  P,j 
-  ET!  !p  M 
E rip -F 

icmtimje 

IF  • I . LE. ?■  THEM 
IP  ' : :.  EC',  n.  .;>1.  oe- 7 n 
v=:fi*:pi  ;• 

c'  E  T  ij  p  f  < 

ELIE  IP  ■  I .  EG.  -i.  OP  .  I .  EO.  E  ■  THEM 


EL  J  I  .  E  '  .  -  •  THFM 

-  t|*T. 

.  r  - .  c  r : 

-  -  j  r.  to 

T  F  •  T  .  a  ‘  ,  tmCM 
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THEN 


N=£ 

GO  TO  14fi 
ELZE  IF  • I . EO. 3 ■ 

N=1 

GO  TO  140 
END  IF 
CONTINUE 
IF  CI.EO.  1>  THEN 
N=3 

GO  TO  140 
ELZE  IF  CI.E0.2.0P.  I.EQ.  3>  THEN 
N=  1 
END  IF 

Z  PH  1  =  1  F'F  ■  N *  KS <  PNX  '.LP-3>  .  PNX  ■:  LP-£>  .  PNX  CLP- 1  >  »  PNX  ■:  LP>  -  X> 
_ PHc  =  H IF1’  I H «  PNX  1  L  1  *  PNX  1  L+ 1  1  >  X  1 
IF  *.ITS.  ED.  c\>  GO  TO  16  0 
IF  CI.EO. l.OP. I.E0.2)  THEN 
V = S  P  H 1  ♦  T  P  H  2  ♦  E  X  P  X  -  H  Z  H  > 

PE TUP N 

ELZE  IF  1  I .  EC'.  ?■'■<  THEN 

V= Z  PH  1 ♦ Z F'H2 

PE TUP N 

END  I F 

CONTINUE 

IF  < I . EO. 1 . OP . I . EC. 6  THEN 

V=  I  PH  1  ♦  I  F'H2 

RETURN 

ELZE  IF  I . EO.  ? '  THEN 
V= Z  PH  1 ♦ IF  H£*> 

END  IF 

P'ETUPN 

END 

FUNCTION  PLF ■ X< I . J> 


♦  THIZ  ROUTINE  Z ELECTZ  hND  CONFUTE Z  THE  hZCOCIhTED  LEGENDRE  * 

♦  FUNCTION  MHICH  II  PEING  INTEGRATED. 


FZ-!  »•:. 


II  w 


IF  ■ I . EO. 0. HND. J. EO. 0 •  THEN 

PLF  =  1  .  0 

RETURN 

ELZE  IF  ' I . EO. 1 . HND.  J. EO. 0 '  THEN 
p  L.  f  ■= : 

C'Etiic  n 

6L  ZZ  IF  I  .  E  ' .  l  .  Hf*  D .  .  .  £  0  .  •  ■  THEN 


CCTI  ;Ui 


z:  :r  i .  “•  . 

*  •  -  +  a  *  p  • 


ELZE  IF  ■  I .  E1*' .  6  .  HN  D .  j .  E  .  1  1  I  HEN 


.  HNp. 
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PLF  =  T'*h 
PETUPM 

el:e  if  '  i. ec.  e.hmd.  j.eo.  cm  them 

P L F  =  1  4 8 .  ♦  1  :•••: +  7 c! .  ♦  X ♦  E  1  4 b . 

petupm 

ELI E  IF  •  I.EO.  E.ftMD.  J.EO.  1  /'  THEM 
F'LF  =  -  1  TOFT  •  ft  •  ♦  >  d yy ♦  X ♦  ♦  c  +  7  c!  ♦  EL1  1  48 
PETUPM 

EL  IE  IF  ■  I.EO.?. HMD. J.EO. £ •  THEM 

F’LF  =15.  *X*Fl 

PETUPM 

EL  I E  IF  1  I.EO.  3.  HMIi.  J.  EG.  1: 1  THEM 
PLF=-  1  5 .  ♦•H**l  .5' 

EMDIF 

PETUPM 

EMD 

I  UP  POUT  I  ME  ftMGLE  1  LMftX .  C  L  <  T  H  «  PL  *  C  ft '  I  IGT*  XTEVL  >  XT  DDL i> I-  T  -  XI  ■ 

D I MEMT  I  CM  C  L  •  0 :  I«  0 :  3  ■  -  I H  •  n :  3 ,  0 :  3?  IE  1  «  PL  •  3*  3*  S'<  «  C  ft  1  3  .  3  >  3  0  1  .  i  !T  E'.'L 
14*.  X .  DEL  ’  U  •  4  1  *  X  .  1  i.i  i  4  1 


♦  THU  PCUTIMF  F I  IEMELE I  THE  TCTHL  hNGLE  IMTEGPftLT  FOP  hLL 

♦  CCMBIMhTICMI  Cc  the  TPHEF'ICftL  hhpmdmic  TF'IftL  ftMD  HEIGHT  fumct: 

♦  IUEXCP  IF'TI  . 


PP I  NT ♦ •  TDThL  IMTEGFftL  VftLUET  I M  MftTPIX  FCPM*  CF  THE  :PhEPIC6L 

C'PIMT*.  hpcmEUu:  TPIHL  ftMD  '"EIGHT  FUN'ITICMI  UTED  IN  THE 

PFIUT*,  GPLEF'  IN  I  GLUT  I  DM  CF  THE  EVEM-F'ftF'  I  TV  EDLTZNPNN  EOUhTIC’- 
PPIMT*.  THEIE  VFLUEI  ftF'E  'ELECTED  PRODUCT'  DF  THE  hTTCCIhTED 
F'PINT ♦  »  LEGENPPE  PGLVMOMIftLT  HMD  ftUHULftP  IMTEGFftL  I  HH I r  H  i.EP'E 

F'F'  I MT  ♦  «  COMPUTED  EhPLIEP.  THEY  PEPPEIEN'T  THE  I C  ft  T  TEF  I  MG  ftMD 

F  P  7  N  t  ♦  «  MCM-  I  C  ftTTEF  I  MG  ftNGUL  ftp  IMTEGFftL  t  EFMELI  WHICH  HPE  FCcmE 
ppi nt*.  uhem  :cHEPicftL  hhfmcmic:  hfe  u:ed  ft-  the  hmgulhf  tpih_ 

PPIMT*.  HMD  '"EIGHT  FUNCTION:  IM  THF  GftLEFT  IN  IDLUTICM. 

F'F  I  MT*. 

F'F'  I  MT  *  . 

DC  EE 0  I  =  I< 4° 

I  l.l=i.  ■ 

rc  •  :x  rri .  i  them 


gc  re  1 1.1 

EL  I E  IF  1  I  .  E 0 .  E  1  THEM 
I  IH=E 
I  F'H=  1 
GO  I C  10 

el : e  if  ■  i . r " .  •  ■  them 


:.c  TlI  1  i.i 
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■50  TO  1  0 

EL  5h  IF  '  I .  EG.  t  .  THEN 

i  :h=5 

IPh=£ 

GO  TO  10 

ELSE  IF  ' I. EG. 7.  THEN 
ISft=3 
I  F'H  =  £ 

GO  TO  1 0 

ELSE  IF  a.  EG.  S''  THEN 
I  S'ft=5 
I  F'ft=E 
GO  TO  SE¬ 
EL  SE  IF  'I.  EC'.'?'  THEN 
ISft=- 

IPft  =  4 
GO  TO  1 0 

ELSE  IF  ■  I .  EG.  1  ’?  1  THEN 

1 3  ft=G 

I F'  ft  =  G 

GO  TO  1 0 

END  IF 

GO  TO  SG 

do  so  i  =o. lmh: 

DO  SO  N=0.1 
I T  =  n 

I  i,l=  I  i,i+ 1 

do  so  l=o.lnp: 

DO  SO  M  =  0  »  L 
IT=IT+1 

f  ft  •  I  M  •  IT.  J'i  =f L  1 1  •  N  1  ♦.*_  L  1  L  .  M  1  ♦  I  H  •  M «  N .  I  _  ft  •  ♦  PL  1  IN*  IT.  I F  ft  1 
C  ONT I  N'.'E 
GO  TO  SS 
DO  SI  f  =0,  lnh: 

DO  35  N  =0W 
I  T=  0 

I  !,)=  1 1,.+  J 


:  t  : r  - 

ft  1  I >  iT  > 
CONTINUE 
IF  ■  I  .  LE.  ? 
IF  1  I  .  F . 
DO  4"  '  -  0 
DO  4 1 1  m=g 


SO  4  i  L -  0 


H  •  N  *  •'i  *  i  ..  H  '  ♦  ^  L. 


i  t  *  :dh  1 


OP  .  I  .  EG.  1'?  •  GO  TO  £  0  f , 
l.GF.  I.EO.E'?'  THEN 
LN4 


"ENF  t  =  1H  ■  !■  •  •  - 

■'?  ■  IT>  I  '  -r L  1 '  •  N  ■*>'!.■  L  . 


*  ,  TH.C  1  *pi  ,  I!,:  ,  IT, 


-n  r1'  ai. 


I .  EG. 
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TEMP  1  =PL  I  l.l .IT-  7>  +PL  I W  -IT-  3> 

CH  I  l.l ,  IT.  D  =CL  CK  .  N>  ♦CL  ■:  L  *  M>  ♦TEMPI  ♦SH  ■  M .  N .  6> 

END  IF 
CONTINUE 
GO  TO  £00 

ELSE  IF  CI.EQ.££>  THEN 
I  SR= 3 
IF'H=3 
GO  TO  45 

ELSE  IF  a.EC'.£3>  THEN 
I  SR= 5 
I  F'H=3 
GO  TO  45 

ELSE  IF  ■ I . EG. £4 .  THEN 

ISh=6 

I  F'H=5 

GO  TO  45 

END  IF 

GO  TO  55 

do  5  0  h  =o. lmh: : 

DO  5  0  N  =  0  •  T 
I T  =  0 

1 1,|=  1 1,|+  i 

DO  5  0  L  =  0.  LMH! ' 

DO  50  M=0.L 
I T- I T+ 1 
M  >L-M 

R  =  1  1  .  -  1  - 1  .  ■  ♦♦MS  1 
I  F  1  M .  E  i. 1 .  0  '  H = h  •  £  . 

CFi  •  I  l>l .  IT.  I  >  =  f  L  '  L  .  M  •  ♦♦£ ♦C L  '4  .  N  1  ♦  S  H  '41 «  M .  I S  H  1  ♦  F  L  1  I M »  IT.  I  F'  R  1  ♦  H  ♦ 

IF1  I . GE . ££. PND .  I . LE . £4 1  GO  TO  £00 

CONTINUE 

IP  '  I.  EC.  £3'.'  THEN 
DO  0  0  I  =  0.  LMH: ; 

DO  HO  N  =  0 .  h 
I  T  =  0 
T  i,i=  I  l,i+l 

SO  t  0  •_  =  ii .  lM FO¬ 
GG  -  0  M-m. L 
I  T  =  I  T  +  1 

m:=l-m 

H  = 1  1 .  +  1  —  1 .  ' ♦♦M :  ■ 

IF  1  M.  EC1.  0  ■  H=P  £. 

CH  1  I  i.i.  IT.  I  ■  =C  L  •  L  *  M  ■  ♦*£♦•:  L  •  V  .  N  ■  ♦  C  H  >  M .  N .  R  •  ♦PL  '  IN.  I  T  -  *.♦*♦'-■.  I  •  L 
GO  TO  £  0 o 


1  OS 


THEN 


I  SB =5 
i:ri=7 
ISE=3 
IP  FI  =3 
IF'B-3 
h 0  TO  130 

else  if  ■  i . E'"> .  i ^  • 

ISFt=3 
1 3  B = E 
I  3D =7 
i:e=s 
IPR=5 

IF  1-3 

GO  TO  130  _ 

el:E  if  a.EU.13'  THEM 

I3Fi=5 
1 3  B  =  3 
I  ED =3 
ISE=7 
IPFt=3 
I  PE' =3 
GO  TO  130 
ELSE  IF  'I 
I  Sh=5 
I  SB =5 
I  SB  =  3 
IS  E  =3 
IFF-? 

IFB  =  3 
60  TO  130 
ELIE  IF  1  I 
I  sp=5 


EC1.  1 4  1  THEM 


EC.  15 


them 


i  :  b=e 
i :  r=3 
i :  E=s 
I  P  Ft  =  5 
TPB—  " 


i:t=3 
i :  B-3 
i  :e=t 

I  F'P  =  S 
I  F  r  =5 


IPB-® 
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I  3£=t 

i  :d=s 
i  :e=8 

IF'H=5 
IP£=5 
GO  TO  130 

EL  IE  IF  ■:  I.  EO.  THEM 
I  IR=F 

i :  e=f 
i :  d=s 
i:e =8 

IF  F  =  F 

IFB=6 

GO  TO  130 

END  I F 

GO  TO  150 

DO  145  1  =  0 «  LMFX 

DO  145  N =  0  > 1 

I T  =  0 

do  145  l=o,lmf:: 

DO  145  M=  0 * L 
I P  =  0 
I T= I T+ 1 

l_  H  1  1 1>1 ,  IT,  I  '  —  U  . 
do  145  jF  =  o,LnR:: 

TENF=0. 

::i  u=: : nr l  '  jf ■  > : igt-":cdl  •  jf  *  > 
if*  i.EC-.E0> u=:::evl*  jf  * 
do  i4n  j:  =  o-jf 

IF-IP+1 

Ift-4 

IF'  J I .  E ' ?  .  0  I  ft  =  I ft  ■  3 

T EMF  =  I  h •*.  L  1  1  » f i  1  L  1  L  '  1:  *  ♦!.  L  *  JF  *  J  .  1  ♦  ♦  c ♦  1  _  H  1  J  _  *  H  *  I  _  ft  1  ♦  _  H  1  j  .  <  fl *  I  _  3  1  — 

i  *  j :  <  tj « i : d  ■  ♦  :h  ■  j :  *  m«  i  :e  ■  ■  *fl  *  it.  if-  iff •  ♦f-l  *  im«  if*  iff  *  +temp 

G  OFT  I HUE 

Cft  '  I  M  *  IT,  I  ■  =TEMF  IJ  +  CF  *  II,',  IT-  I  1 


i  :f=i 
i  :e=g 
i  :  r  -T 


i  :e=': 


r  i  "  f  ir  *  cr.  t  j  ,  tlch 


i  :e=  - 


ELIE  IF  •'  I  .  EC1.  £7>  THEN 

I  1H=F. 

I  IE-F. 

I  ID =8 
I  .  E=F; 

IFR=5 
I  F'E-6 
ENDIF 

DO  175  r  =o<  lnh:  : 

DO  175  N=0-K 
I  T  =  U 
I  M= I M+ 1 

DO  175  L=0-LMHK 
DO  175  M=  U< L 
I F'  =  0 
I  T  =  I  T+ 1 


C  R  1  I  FI »  I  T  «  I  1  —  U . 

do  175  jf=m.lnr:  : 

TEMP=U. 

: t  u=: : : odl  •  jp  ■  ■  igt-:::odl  ■  jp >  -• 

DO  170  j:  =  o.jp 
IF-IF+1 

m:=l-m 

R  =  E*  ■  1  .  -  ■  -  1  .  •♦M  I  ■ 

I F  -  J  I  .  EC' .  0 .  OF  .  M .  EC  .  0  >  H=H  £ 

TEUFEL  ■  L-  M  ♦CL  1  h  -  N  -  ♦CL  1  JP-  JI  •  ♦♦£♦•:  IH  '  J  C  .  H-  ICR  ■  ♦  EH  '  J  I  .  M  <  I  IE 

i-ri.  i : d *  ♦ : h •  j: .  r-u  r  :e'<  ■  ♦fl  < it- ip.  ipr.-  ♦pl  -  iw.  ip.  ip?  ■  ♦r+temp 

CONTINUE 

CR  •  I  '.I  -  IT  -  I’  =CR  •  I  l.i  -  IT  -  I,'  +  TEMPVT  IJ*KC  •  L.1 
CONTI  MCE 

PRINTS-  NRTF'  I ' '  VRLUEZ  FOP  RNGLE  INTEGF'RL  =  -  I 
DO  El  0  Il.i=  l  ,  f  T 

PRINTS-  '  -  CR  ■  I  l.i  -  IT  -  I  •  -  I  T=  1  -1  J  • 

PRINTS. 

F'F  I  NT  *  - 

CONTINUE 

RETURN 


♦  THi:  F  CUT  INF.  COMPUTE:  THE  EPWEPICRL  HRRt'ONIC  EMPRNS  ION 

♦  coefficient:  op  ndfmrlizrtion  conctr-t:. 

PERL  INF. IDF 


4 


;n  th 


I  DF  =  I DF ♦ I 
FT= INF  IDF 
GO  TO  100 
FT  =  1 . 

CL ■ L- M '• =  7  OPT ■  ■ E  ♦  L  +  1  ■  ■ 4 ♦ ■ 4*mThn >]...■ t 

continue 

RETURN 

END 

S  UFR  CUT  I  HE  ."PL  I  HE  ■  7  P! 1  -  N.'  7  <  TDL  •  f  OUNT  ■ 

D I  MENTION  PNX  •:  :  7  '•  *  C  PL:  <  7  •  7 «  E  • 

COMMON  Hi  FTP' >  Ih  3  JF  -  JT  F4  L  » 1 


♦  THI_  ROUTINE  COMPUTE^  THE  FICUFIC  IFLINE  INTEGRAL  7  BY  CALLING  ♦  ♦ 

♦  ft  TEN  POINT  NEWTON-COTE!  INTEGRATION  ROUTINE. 


DO  9  u  JF'=  1  i  r PL 
DO  90  JC  =  1  -  N7  YE 
-  F’  !■■■! J  F'  ’  _l  i_  »  I1  —  m  .  ij 
JT=  JC-.JR 

IF  '  R  pi  1  J  T  1  ,  i  b  T  .  1  bD  TO  9  0 

I  F  =  4 

IF  ■:  JT.  EC.  O'  THEN 
F  1  =  1 
GO  TO  5 

EL  7  E  IF  <  JT  .  EO  .  1  ■'  THEN 
F  1=3 
'30  TO  5 

EL  7 E  IF  ■  JT.EO.E'  THEN 
1  1  =  7 
GO  TO  0 

EL  7  E  IF  1  JT  .  EC1 .  3  '■<  THEN 
F  1=4 
GO  TO  5 
END  IF 
F  1  =  1 

IF  '.'JT.EO.-l  ■  THEN 
F  F  =  3 


GO  TO  f. 

EL  7  E  IF  'JT.EU.-3'  THEN 
1  F=  1 
GO  TO  0 
ENDIF 


conttni  :e 


i  •> 

i  _ 


1'  1 

IF  -JP.E0.4'  THEN 
IF  'IP.EO.l'  THEM 
IFF  =  Z 
FETMPM 

EL  ZE  IF  <IP. EC. E  '  THEM 

:pf=zi 

FETMFH 

ELIE  IF  *  IF. EM.  ? 1  THEM 

:ff=z+:  >zi 

PETL’PH 
ELD  IF 
EMT  I F 

*i  =  Z  — -t ♦  F  ♦♦  ? 

V  I  =  t  Z  .  ♦E:*-*£  +  Zl 
IF  'JP.E0.3>  THEM 
IF  ■ IP.EO.l >  THEM 
IPF  =  Y 
PETMPM 

el:e  if  'If.eo.e,.  them 

:  pf  —  'i  ■  i 

petijfm 

EL  "E  IF  ■ IP.Fm.?-  THEM 

■PP  =  ',  4-  1 


em:  :f 

V  =  Y  +  P  *i"  *♦  T 

t  •'?  .  ♦  »!  ♦  ♦£’ 

rp  •  _>F. EM.  £  ■  THE 
I F  •IP.EO.l'  THEM 


c  r  t 


in 


i'=: :  l 

l  .1 = V — 4  ♦  ♦  ♦  5 
Ml  =V1+1 

IF  • JP.EO. 1*  THEM 
IF  • IP. EM. 1'  THEM 
ZPF=l.l 
FETUPH 

el:e  if  •ip.eo.e-  then 
:  f  f=i.i  i 

PETUF  H 

ELZE  IF  ■  IP.  EO.  3  ■  THEM 

:pf=m+:  >i.ii 

EtHD  I  F 
EHDIF 
FETUPH 
EHD 

ZUBPDUTIHE  lOUPCE  1  I  VM«  HI  IT  •  TOL<  FOUNT  • 

dimehz  ion  ivy  •  r-  ?.*  z  -  •  p re:  •  -a:  r  • 

CDMMDM- Hi  FfP. I  H4  1  Z. IH-LF  FJ  L-l 


♦  THIZ  POUT  I  ME  ZEIECTZ  HMD  COMPUTE  I  THE  IHTEFPOLhTED  ZDL'PCE  * 

♦  IHTEGPhL  EV  OH.  LI  hi?  H  TEH  FDIHT  HEUTDH  COTEI  IHTEGFfiT  I  OH  FOOT  I  HE 


do  so  lp=i»hyi 

DO  3  0  L C  =  1  «  H'  I-E 


VI  !  1  L  F  « l_iw  1  I  1  ~  1 1 .  i.' 

LY=LC-LP 

IF'LY.GT.  1.  OP.LY.LT.  -  Z:  •  GO  TO  30 
1  F  =  4 

IF  'LY.EO.O  THEM 
1  1=4 

GO  TO  1 0 


ELZE  IF-LY.EO.m*  THEM 
1  1=  ■ 


GO  TO  1  0 
EHDIF 

IF  1  L'i  .  EC1 .  - 1  •  THEM 


ELZE  IF  •Li. EM. -3'  THEM 
1  1=1 
1  F=3 
GO  TO  in 

EL  _  E  IF  'Li.  E'"1.-:'  THEiM 


tu=  i 

I  0  !  ■  ■  -  '  -r  :  .  '  F 

~ ,  r;  _  ~ 

L  E' =  L  F  ~  M !  .  f  . 

V*!.. ! .  '.;h  i  '  M  '  *  L  i  <*  •  1 .  -  '  «  '  ;..  .  *  *  '  r<  1 
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:  vx  'LR.  lc  •  n  =svx  •  lf. lc  .  i  >  +tv 

I H= I H+ 1 
CONT I  NUE 

PRINTS.'  M  PI  TP  IK  SOURCE  INTEGRAL^  <  I 

do  ?o  mp=  i  .  nxs 

PR  I  NT ♦  <  1  '  ■'  »  i'VK  1  UR .  MC »  I  >  *  HC=1  *  NXS  — 2  > 

PRINT*' 

RETURN 

END 

FUNC  T I ON  H I F -  I H - K 1 ' X2 » K> 


♦  THIS  ROUTINE  COMPUTES  THE  LINEAR  LAGRANGE  POLYNOMIALS  WHICH 

♦  ARE  BEING  INTEGRATED  AS  A  PART  OF  THE  SOURCE  INTEGRALS. 

IF  1  I H . EG. 2>  THEN 
H  I  F  =  K2-: :  <K£-K1> 

RETURN 

ELSE  IF  •  IH.  EC.  1  :■  THEN 

h  i  f=  :  • 

END  IF 

RETURN 

END 

_  IJ  B  F'  D  U  T  I N  E  DIRECT  1  L  H  A !  N  A  F  D .  FHZ'  PNP '  N  P  S  »  NO.  >  S  I G  T  •  HB  »  Y  D*  ASH  • 

D I MENS  ION  ARD  ■  5  •  5 »  S  '  .  PNZ  ■ -2 :  ?"<  .  PNP  1  -2 :  7  < 


♦  THIS  ROUTINE  COMPUTES  THE  DIRECT  SOURCE  POINT  VALUES  WHICH  ARE  * 

♦  USED  IN  THE  FIRST  COLLISION  SOURCE  INTERPOLATION. 


PRINT*.' 

PRINT*.  INTERPOLATION  POINT  VALUES  .  ELEMENT  <  I  <  J'> 

1  .  R' 

PRINT*. 

11,1=0 

do  so  l= o . lma: : 

DO  5  0  M= 0 .  L 
11,1=11,1+1 

DC -40  IT=i.NZS-S 
:c  4c  :=i,NF:-s 
ZH— 'MZ  -  I  r -  i  , 

r’^H  =  .  vF  T  •  r  NR'  •  I  —  [  1  •  *.-+_H**c 

IF'F'^H.EL'.  1 1 , !  t  O  TO  A  U 
UD  =  ZH  R'ZH 

IF'UD.EO. O'GO  TO  20 

T  H  l  J  =  _  I  U  T  *  •  E  X  F'  1  —  H  B  A  .  H  1  —  E !  !  F'  1  —  F'  N  Z  1  I  T  —  1  1  A  .  H  1  '  *  A  .  H  U  D 
GO  TO  SO 

TR*l  '=  ~  I  iS  *G'  'C  ,  -c-fj-1  ,  ;t_,  .  ci'u  .  *cfic-  ,  y  —  i  . 

;  r  ■  :  •  p*E.-  ’F  1  -T  AIJ  .  •  '  <-  .  » ,->  r u  •  M  .  ■  .  +  4  —  ■  *  -  !  P 

L  f  •* _ _  •  I  T  —  L  -1  .  H  ■ 

'"OUT  [Mi  ip 

••  4  t  r 1  t ■» .  '  "h  tr  ; ;  '  _’E  "  a n 1 1  r R -  r  f_'i  7  ”  Z'1  ~  'S  I !"  T  _ 


■-  F  I  * :  I  ♦  ■ 

DO  4^  MF  —  1  - v ' Z  7 -2 
-  c-  r  f  i  r  * .  .  .  r- 


FOP  PD  I  NT 


MU  ,  . 


b D  TO  7 U 

PRINT*.  IMPROPER  INPUT  DhTh.H  NODE  CANNOT  EE  AT  THE  EUR  I T  FOINT-Z 
1=  •  PNZaT-l' 

:  TOF 

RETURN 

END 

SUBROUT  I NE  E'VEC  ■'  PiF'D .  SVZ  -  3  VP .  IE.  NZZ .  NR  S' .  LMAK .  CA  .SH3T.PT.  NT'* 

D I M  E  N  _  ION  A  P'  D  1  S .  7  .  S'  .  .  V  Z  1  7 .  7  .  ■■ .  •  .  _  V  F1  1  7 . 7  .  .  _  E;  1  7  ij  1  .  i_  A  *  3 .  i.i  > 


*  THIS  ROUTINE  ASSEMBLES  THE  PPOELEM  SOURCE  VECTOR. 


PRINT*. 

PRINT*.  SOURCE  MATRIX  COLUMN  VECTOR.' 
PRINT*. 

13=0 

DO  9  0  .JZ  =  1  .  NZS 
D  0  9  U  J  P  =  1  .  N  F’  i 
JH=  0 

DO  SO  JL=0.  LMh:  : 

DO  SO  JM=i).  JL 
JA  =  JA+ 1 
I S  =  IS  *  1 
SB  ■'  IS  '  =0. 

DO  77  N  = 1 . 7 
S  C  =  2  .  0 
TEMP  3=0. 

I A  =  i ) 

DO  7  0  IL=0.LMH': 

DO  70  I M= 0. IL 
I  A=  I  A+  1 
TEMPS- 0. 

DO  SO  IZ=1.NZS-S 
do  so  ip=i.np:-s 

IS  R  =  I  R  —  j R 

IF  '  I L  Z. AT .  1  .  OR  .  I P  Z .  LT .  -  S'*  GO  TO  SO 
IF'  TOP.  ST.  l.DR.  KF.LT.-3 -GO  TO  SO 


N^r  37 
90  TO  N'i 

el:e  ip  -to pc. s-  then 
NZ=  ! 

.up  = 

T  T  *  «  .  '  ^UCC 


L 1 0 


r  cr 


■  r 1 .  p  .  i  ■ 


TL.crn 


GO  TO  50 

ELSE  IF  •  N. EG. 5 ■  THEM 
HZ=  1 
MR =2 
Oh  =  26 
GO  TO  50 

ELSE  IF  • M. EG. 6>  THEM 

M2 =2 

MR  =  3 

MR=27 

GO  TO  50 

ELSE  IF  a1.  EG.  7  •  THEM 

NZ=3 

MR  =3 

MR  =  23 

END  I F 

TEMP  I  =  SVZ CJZ .  12.  f  IZ.<  ♦  i-  VP  '  JF  .  I R « f  #F'.>  ♦RF'P  I Z  •  IP  •  I  R,*  ♦C  R  1  JR  .  I R «  MR.-1 

TEMR2  =  TEMR  l+TEMF'2 

COMTIMIJE 

TEMP3  =  TEMF'2+TEMP3 

IF '  MR. EC. 23. OP . MR. EC. 26>  S C  = - 1 .  0 

IF  CM.  EC.  7">  GO  TO  72 

T  E  M  P  3  =  T  E  M  P 3 ♦ 3 ♦ R  T  R M -  1 .  <  ■•  3 IGT*SC 

GO  TO  75 

TEMP  3 =  T  E  M  P  3  ♦  3  ♦  R  T  R  N -  1 . > ♦IC 
SB  ■  IS.'  =TEMP 3  +  3 B  •  I3> 

COMTIMIJE 

F'RIMT  ♦  «  « f  T  i  ELEMENT  '-O  TO  ELEMEMT  NUMBER*5 "  >  I 3 

PR  I M T ♦ ,  <-  ,  S  B  <  1-  S  >  .  h  !  =  i  :  -  '  h  T  -  1  >  .  I  S  ' 

COMTIMIJE 
PR  I MT ♦ . 

F'RIMT*.-- 

PRINT*.' 

RETURN 

END 

SUBROUTINE  STIFF  ■  CR.  3  PR.  SFZ.MZS  .  MR S  .  S  I G T .  R S -  PMR .  F'MZ -  LMR!  ' .  NT"' 

D  I  MEM':  I OM  L  R  '  U .  r’ .  0  ’  .  ;.  f~‘P  •  7 . 7 . 6- '  *  _  F’Z  1  7 .  7 . 5  1  .  R  _  ■  5  0 .  5  0.'  ■  F'MP'  •  — c  •  •  1  .  R 
MZ  1  -  i  '■ 


THIS  POL’ I ME  R'SEMBLEI  THE  F PC EL  EM  COEFFICIENT  OR  STIFFMEII  MRTF I 

BO  2  0  J  = i *  M T 
BO  20  f =1 . MT 

r  r  ■  j  •  i  '  =  o .  o 

I  J  =  0 

BO  ' Z  ’  •  •  •  '  . 

:  0  -  I  7  =  !  .pf 

:  ,  f  f  ■  |  Ul  ’  * 


11 


ELSE  IF  CN.EO.ll>  THEN 
NF-£' 

HZ=  1 
NE=1 

GO  TC  30 

ELSE  IF  CN.EO. 1£>  THEN 

NF-4 

NZ=£ 

go  to  so 

ELSE  IF  CM.  EC1.  ISO  THEN 
NP=£ 

NZ=  1 
NE=  1 

GO  TO  35 

ELSE  IF  CM.E0.14'  THEN 
NR  =  3 
NZ~  1 

GO  TO  30 

ELSE  IF  cN. EO. 15>  THEN 

NR  =  5 

NZ=£ 

NE=1 

GO  TO  30 

ELSE  IF  'N.F.O.lt'  THEN 

NF-4 

NZ=3 

GO  TO  cl 

ELSE  IF  cN.EO. 17>  THEN 

NF-5 

NZ  =  £ 

NE=1 

GO  TO  33 

ELSE  IF  CN.EO. 13>  THEN 
NF  =  6 

mz=  s 

GO  TO  30 


Itf.  ]  G  •  Tl-r*  1  '  .  ■  *  -  I'bT 

i  H  «  .l  '  *■-  *H  T  HN  '  1  .  -  i  tNF'b. 


-  S-  ,  -  /.  r  .  '  .  f3c  .  I S  .  -  T 

■  v  .  ~.r  ■  . ;  'P  ■  .  GT .  1  '•'0 
-r  /  1  ■  THEN 

1.1  —  4 


riu  li 


.p  ,  NP  '  *  I  r  S  '  I  Z  ■ 


I-1.. 


'CM,  -.CO.  1 5b  NR'*  ♦SFZ'IZ*  JZ  •  M_-' 

GO  TO  S3  _  _ 

TEN=:FR  •  IP.  JR*  NP">  ♦SP_  ••  J-  ’  I  —  ' ' '  —  ‘ 
i.D  TO  ;:c!  _  T  _  . ,  — 

TFN=-pp'  IP.  IF.  NR"-  ♦IFZ-  J-  ' 
t|!,F  j  ,TEM>CH  -  JH.  Ift.  N  ■  ♦  ■  "I  •  ♦♦NE+TEMPl 


■CMP|  = TEMPI  •  -  ♦H  :  HN  •  1 

"  P'  M  C  S  ~  ‘  P’C  I  P  «  [  ^  H  1  ♦  . 
--Mr  ~  Cl  ,  •  C-  ,  TC,r  «  ♦  ’ 


; .  a  .  ♦>  -t  1 
■ ,  4  >  ♦ H  • 


m 


el:e  if  ■  Li'. eh. -i ' 
L  E  =  7 

I  rrji 

,-Q  T;  C" 

EL  IE  IF  'LD.EE.-Z 

.  THEM 

LE  =  E 

ED  TC  75 

END  IF 

ELZE  Ic  '-'ZZ.EC'.E' 

l  i  = : 

IF  >  Lr.EO.O  '  THEM 

THEN 

EC  TC  7  5 

EL  IE  IF  'LI'.EO.l' 

THEM 

LE  =  4 

ED  TD  75 

ELZE  IF  ■ LD. EE. ~ 1 
LE  =  E 

■  THEN 

EC  TD  75 

EMI' IF 

ELZE  IF  •  jZ.EE.  Z- 
L  1  =  E 

IF  'Ll  .  EE.  O'  THEM 

then 

Li  -E 

EC  TC  Z  5 

ELZE  IF  'LH.EE.!' 

t  U  E  N 

LE=  7 

EQ  TD  75 

ELZE  TF  'LH.EE.E- 
LE  =  4 

END  IF 

T  HE '  ■ 

EfJD  I F 

FEZ  Z  =  Z F f '  1 . L  1  *  FHZ 

■ ,  j2-  • .  =  nz  ■  j; 

I LE ’ PNC  1 IZ-7  ' *  PNC  1 

I z  — E  ■  '  fnz  •  IZ- 

ED  TC  4 11 

!■  i‘=rc  z  -  jz 

►■I  T.  -f  J^Z  -  I  Z 

TC.i  r, .  |-r  .  PP.  ME. 

i;t.  ,  nn  Tn  j. 

.  t-  •  • 

- "  . 

’  ,  ,  -  r  r  - 

;■  :  Z  -  ’ ;  -  ’ 

'  z 

»  D=MF  Z-JP 

r.'E-r.c  ~  -  ip 

’  ‘  H.  'T.  E.  DF  .  MI' 

.  i?T  .  c  •  bU  TG  4 

I  UBF'OUT  i  me  eg1- HD  ■  j'  :•  i: fn;  ,  temp.  m: :  :  ■ 
D I  MEN  -  icu  fm:  :  •  r  - 


♦  thi:  f out i  me  '.elect:  hud  compute:  the  vhlue:  of  the  bicubic 

♦  :  PL  I  ME  FUUrTIOMZ  HT  THE  OUTEP  PHD  I  US  HMD  TOP  ZUPFRCE  OF  THE 

♦  F  POE  LEM  C'CL  I MDEP .  THE  IE  VRLUET  RPE  PEOUIFED  I M  THE  EVRLUHTIDU  OF 

♦  THE  .  ijF  FHl  E  IMTEhPRL  TEF’Mo. 


l  z=r<: : :  -  j: : 
ld=  j:  I 

IF  'Ll. EC. O'  THEM 
L  1  =  1 

IF  •LD.EC.O'  THEM 
Lc'=  1 

CO  TO  4n 

EL  I E  IF  •LD.EO.l'*  them 
6 0  TO  40 

el:e  if  •ld.ec.E'  them 
le=  : 

60  TO  40 
EMDIF 

EL  IE  if  'LI. EC.  I  '  THEM 
L  1  =B 

IF  ■  Lr.  EC'.  0  •  THEM 


60  TO  40 

ELZE  IF'LD.EC.-! 1 

THEM 

LE=  1 

60  TO  40 

ELZE  IF  ■ LB. EC .  1  ' 

THEM 

le= 

60  TO  40 

EMDIF 

ELZE  IF  'LZ.EC.E' 

THEM 

L  1  ~  0 

IF  ' LD.EC.O'  THEM 

L  E  =.£ 

60  TO  40 

ELZE  IF  •LD.EC.-E'  THEM 
LE-! 

EMDIC 

ChT|  T  C 

rc-MC  -  -fc  ■  i  .  I-':-’-  •  :  ■  "  ‘  -  -  ,  .  c  o  ,  ■■  -  <  ,  n."  ,  .  .FM 

-  ■  ■'  "  -  ••  >  .  r  :  '• 

F  F.  T'_IF  M 


PRINT*.  •  MR  TP  I Y  CHECK  FOP  DIRGONRL  DOMINRNCE. 
PRINT*. 

PRINT*.  ROM  D I RGONRL  ELEMENT  CUM  OF 

PRINT*. 

DO  40  1=1. NT 
TEMF-0.  0 
DO  30  J=l. NT 

temr=temp+rbc  (hs  ( i .  j:>  • 

TEMF'  =  TEMP-REC  '  Hi  I  .  I 

PRINT*.''  .  I  .  ■  .  H  j  1  I  .  I  1  .  .TEMP’ 

PRINT*. 

CONTINUE 

PRINT*. 

PRINT*.  MR TR iv  CHE  Cl  FOP  CVMMETPY. 

PRINT*. 

DO  CO  1=1. NT 
DO  CO  J= 1 . I 

F'P  I  NT  * .  H_'  .I*  .  'J.  •  =  .  H  _  *  I .  J . 1  . 

1  C  (  J.  I  - 
CONTINUE 
RETURN 
END 

CUEROUTINE  I OL VE • N . 3 E . R ‘ . NT . !  FT 
D I  MEN  3  I  Of4  R :  •  5  0 . 5  0  ■  .  3  E  ■  0  0  1  .  :  H  •  C  0  ■  . :  E  ■  5  0  • 


CUM  OF  PON  element:. 


.  J  .  .  .  I  . 


*  THII  ROUTINE  ICLVEC  THE  PROBLEM  MRTR I C EC  BY  THE  ITEFRTI"E  MET 

*  OF  ZUCCEICIVE  ovep-relr:- rticn. 

FPINT*. 

PRINT*,- 
DO  .3  0  1  =  1.  NT 

:  :r - 1  ■  =o. 

:  :e  ■  i  ■  =o. 

I  T  =  0 

I T= I T  + 1 

DO  j  0  1=1, NT 

temp=:e  •  i  •- 
”  ’=;.nt 

•rvDsTCWp-C,:  .  1  . 

-  •  ’•  -  rCMC  »ii!  ■  P:  T  '  +  "l  1  f  1 

DO  CO  1=1, NT 

IF  1  'E  1  I  '  .  EO1.  0.  >  CO  TO  C  0 

TOL=  '  '  P  '  I  ■  -YR  '  I  '  '  -YE  •  I  • 

IF  1  HE  :  '  TGL  1  .  GT.  .  Ii01  '  CO  TO  7  0 
continue 

i-O  tv  or, 


DO  in  I  =  1 . NT 


“  c  I  *  *  •  r;jc  “  CL  N  T  T  Cl-  n  H  :  :*  !•..  ere  COD  .  THE  ^  I f 

ff;nT*.  the  number  of  iteration:  MERE  .IT 

lO.urrr!'  ••FC  — .  C;  PH"  "  ION  COEFFIC  IF 


PRINTS.  THE  PPOELEM  HRG  NOT  CONVERGED  IN  -IT-  ITEPHTIONT. 
PRINT*. 

PRINT*.  THE  ERROR  =  .  TDL 
PRINT*. 

PRINT*.  >  . XE • I • • 1=1 .NT ■ 

RETURN 

END 

I  IJEF'OUT  INE  TFLUX  1  PNR .  F'NZ  .  XB .  R .  Z>  1  T.NZS .  NRS  .  TPHI .  DPHI  .  HE'.  I  I G  T  .  R  C  H  • 
ID' 

DIMENSION  PTC  ■  7"<  .  PNR  •  -£:7‘>  .  XE  •:  0  0  >  .  I  CPF  ■  4-  4 

*  THU  ROUTINE  COMPUTES  THE  GCRTTERED  RND  TOTRL  REhCTICN  RRTE 

♦  fluence. 

CRLL  IEMftT  •.  I C F'F .  PNR «  F'NZ .  F  .  Z-NRS  .  NZI  :■ 

1  Z  =  5 
TEMP  1  =  0. 

DO  40  1=1.4 
MZ=I  I  F'F  1  I  .  1  ■ 

IF • HZ. GT. NZI ' GO  TO  40 
1  Z  =  1  Z-l 
>  F'  =  c 

Du  I  0  J = 1 . 4 
up= i :rp ■ j. £  > 

IF ■ NR. GT. NRI ' GO  TO  CO 

NF=  1  *>■  T*  '  ■  NR-  1  ■  +NR  I  ♦  1  NZ-  1  '  ■ 

f  R=1  F-l 

TEMF  =:  l  •  NF  •  ♦  F'F  ■  1  . 1  Z.  F'NZ  <  HZ-  C1  .  F'NZ  •  NZ-cV  .  FTC  •  HZ-  1  '  .  F'NZ  '  HZ  '  .  Z  ■  ♦  I  '  r 
1  1.1  F  .  F'NF  ■  NF  -  '  .  F  N F'  •  1*1  F'  —  c •  .  F' N F  1  N F'  —  1  *  .  F  N  F  '  M  F'  1  .  F  ' 

TEf'-Pl=TENF'l  +TENR 

CONTINUE 

CONTINUE 

H=  I  CRT  1  1  E  .  *RTRN  '  1  .  '■  ■ 

TF'HI  =h*TEmP1 
Z  H = Z  -  H  £' 

F  ZH= I  OF T ' F  **C  +  ZH**c ' 

I  F  1  R'ZH .  E ! *  i  '  GO  TO  V  0 


-:-l  ;  =  f  :  -  T  »  ■  P  '  ;c  •  p  <4  *H  •  'P  •  -Z  H  H  «  <  *H  *  H  1  ;  D 
GO  TG  ~0 

T HU =IIGT*E:  F  '-Z  h7H'*R 

DF'H  I  =VD*E:  F'  ’  -TRU  '  ■  1P*RTRN .  i  .  .  *RZH**£  ' 

GO  TO  GO 

PF  I  NT*.  ERROR.  THE  Fi_C  CfiNNOT  EE  CONFUTED  R  T  THE  EUR  -  T  POINT. 

■  THD 

C  - '  K  *  1 

i  ■  .  r, 

,  •peri  IT  INC  tCMu}  ,  *  -  p  r  ,  Cue  .  pN~  .  p  ,  -  .  f.c  •  ■ 


Nl  I  P  :  OF  -~E  r  1CUF  ;c 


*  .ELE1  TED 


pc int : . 


LIU  4 I  I  I  =  H«r)_ 

IF  '  Z.  GE .  F'flZ  '  I >  .  HMD . Z . LT . F  MI 1  I  +  1 1  ■ GO  TD  5  0 
CUNT  I HUE 
DU  60  J= 1 ' 4 
I  IFF  ■  .! »  1'  =  1  + J 
DO  7  0  1  =  0' HP  1-3 

IF-P.GE.  FNP  ■  I >  .  FIND .  P .  LT .  F'NP  I  + 1  >  CD  TO  8 0 
CONTINUE 
DO  S5  J  = 1 «  4 
I  CPF  ■  J'  £>  =  I+-J 
GO  TO  100 

PP I  NT*  ’  THE  CF'HCIHL  POCITION  OF  .  O'  '  OP  -P' 
1PPOELEN  D I MEN  : I  one . 

Z  TOP 
FETUF'N 


IC  MOT  lil  I  TH I  N  THE 


Appendix  K 


Numerical  Results 

The  computer  subroutines  which  are  listed  in  Appendix  J 
were  used  to  produce  numerical  results  for  varying  problem 
spatial  mesh  sizes  and  degrees  of  the  spherical  harmonic 
trial  function  expansion.  Some  of  these  results  are  presented 
in  Figures  10  through  44.  They  are  valid  for  the  problem  para¬ 
meters  which  were  presented  in  Chapter  V.  However,  the  cross- 
sections  which  were  used  do  not  accurately  represent  the 
values  for  air  at  sea  level.  Also,  the  air-ground  interface 
was  not  included  in  the  problem  domain  and  therefore  all 
ground  effects  were  ignored. 

Because  of  the  time  constraints  on  this  research  project 
neither  an  evaluation  of  the  accuracy  of  these  results  nor 
a  comparison  to  a  discrete  ordinate  or  Monte  Carlo  calcula¬ 
tion  was  accomplished.  Therefore,  the  results  are  presented 
solely  in  an  attempt  to  show  that  finite  element  space- 
angle  synthesis  is  a  viable  solution  technique  for  solving 
the  two-dimensional  steady  state  anisotropic  Boltzmann 
equation.  They  are  not  meant  to  represent  a  precise  and 
exact  solution  to  the  air-over-ground  problem,  but  rather 
to  demonstrate  that  FF.SAS  may  be  a  feasible  alternate  solution 
technique  to  Monte  Carlo  and  discrete  ordinates. 
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